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Abstract. The equations governing the quasi-one-dimensional steady flow of a 
conducting perfect gas in crossed, transverse electric and magnetic fields are treated 
under the assumptions that the electric conductivity ‘is a scalar, that the wall drag is 
small and that the magnetic field due to the currents in the gas is negligible. The equa- 
tions are normalized and different flow situations are illustrated by means of phase 
diagrams (drawn for a gas of constant specific heats, y = 1.33). The possibility of a 
smooth transition from supersonic to subsonic motion is pointed out and the phenomena 
of standing shock fronts and choking are surveyed for constant-area, small-friction, 
power-yielding flow. The appendix contains some general remarks on the system of non- 
linear differential equations: 2’ = P(x, y)/R(x, y); y’ = Q(x, y)/R(a, y). 

Basic equations. We postulate steady flow of an (almost) inviscid fluid in a structure 
somewhat like the one outlined in Fig. 1. We assume that we have local thermodynamic 
equilibrium everywhere, that the fluid is a perfect gas, which may or may not have 
constant specific heats in the temperature range considered, that there exists a scalar 
electric conductivity and that heat flow and radiation, gravitation and boundary effects, 
except wall drag, may be neglected. We assume that wall friction is small and can be 
taken into account by employing the conventional gas-dynamic drag coefficient for 











Fic. 1. Outline of flow structure. The conducting gas flows in the channel between the two electrodes and 

two insulating side walls (not shown). The gas velocity, u, and the electric field between the electrodes, 

E, are perpendicular to each other and to the magnetic field, B. If u > E/B, an induced current, J, 

will flow through the external load, R.; if u < E/B, R, should be replaced by a supply voltage, V, and 

the current should flow in the opposite direction, producing an accelerating volume force, 7B, on the 
gas. 
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one-dimensional flow. We also assume that the magnetic field due to the currents in 
the gas is negligible and that, averaging over the cross section of the duct, we may 
neglect second-order terms in the variations from the mean values. 

With positive directions according to Fig. 1, we then obtain the following three 
groups of equations 


a) Gas relations 


The equation of state p/p = RT 

The specific heat c, = dh/dT = Rp, B = B(T) 
For 8 = const. the enthalpy is h = Bp/p 

and the entropy is s = R In (h’ Pp) 

The Mach number M = |u| (Bp/(B — 1)p)*” 

b) Electrical relations 2. 

The total current [ = [ (Az)i dx 

Ohm’s law, external circuit V =E Ay = IR, 

Ohm’s law for the fluid i = o(uB — E) 

c) Mechanical relations 

Conservation of mass puS = const. 

Conservation of momentum puS du/dx = —iBS — Sdp/dx — F 
Conservation of energy puS d(h + u’/2)/dx = —iES 


where iS is the electric power yield per unit length. Across a discontinuity (shock 


front) we must have 


( pu = const. 


| p + pu = const. (1) 
lh + u’/2 = const. 

and, naturally, the entropy of the fluid may not decrease in passing the discontinuity. 
Nomenclature 


MKSA-units are used and the symbols have the following meanings 


Az] 

p pressure Ayr dimensions of apparatus 
Az 

p density I external current 

R gas constant (average) V_ external voltage 

T absolute temperature R, external resistance 

h enthalpy per unit mass t current density 

s entropy per unit mass a __ electric conductivity 

8 normalized specific heat 3 magnetic field 

at constant pressure 7 electric field 
u gas velocity S cross-section area 


M Mach number F wall drag per unit length 
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Principally the same equations have been treated by various authors using different 
simplifying assumptions to render them analytically tractable. Thus Way [7] has an- 
alysed, with numerical examples for power generators, three cases of frictionless constant- 
conductivity flow in constant electric and magnetic fields. These assumptions, obviously, 
do not yield a fully determined problem; thus Way introduces one of three restrictions: 
constant velocity, constant area or constant pressure. Neuringer [3] has analysed the 
constant-area, constant-conductivity flow in constant electric and magnetic fields. 
Especially, he has derived equations giving the electric field corresponding to maximum 
power for given length, magnetic field and inlet values. He presents numerical solutions 
for a certain inlet Mach number, which qualitatively agree with what should be expected 
according to our phase diagrams. However, the asymptotic solution, corresponding to 
u ~~ E/B at the inlet, should not give a finite limiting length (6,,,...). Berner and Camac 
[1], when considering a space propulsion device, have studied frictionless, constant- 
conductivity flow in a channel where the cross section and magnetic field are constant 
(but where the electric field varies). As further restriction they choose either constant 
temperature or a certain constant ratio between applied and induced electric field 


(E/uB (8 + 1)/2), which they somewhat arbitrarily designate as the maximum 
acceleration case. This case has also been analysed by Resler and Sears [5]. Sutton [6] 
has made a study of frictionless flow in a constant magnetic field and with E/uB = 4 

: g 2 


(corresponding to a local matching of external and internal load, i.e. iE = 7°/c). Gener- 
ally, he assumes constant conductivity, and as further restriction he uses constant 
velocity (also for varying o), constant temperature, constant pressure, constant density 
or a cross section of constant area but varying shape. The slightly different problem 
where the currents close in the gas has been analysed in connexion with two experi- 
ments by Patrick and Brogan [4]. Finally the constant-area flow with friction in constant 
electric and magnetic fields was studied by the present author in a preliminary report [2]. 
In all these reports the gas was taken to have constant specific heats. 

Normalization of the problem. In the following we shall treat the problem in a 
certain normalized form, which will enable us to give a qualitative outline of the possible 
types of flow by means of phase diagrams. The normalized problem may also be prefer- 
able for numerical integration. 

We denote by a and b two characteristic quantities of dimension length and velocity 


respectively 
a = pu/B’ ofp 
b = E/B 


° ° ° 1 
and let subscript 0 indicate some chosen constant values . 
We introduce two dimensionless variable: 


7) 


u* = u/Ddo 


% 
lI 


Bop/ pub 


and two dimensionless parameters: the channel divergence parameter 
‘The quantity 8) appearing in the equations is unessential and may well be put equal to unity. 
To obtain conformity with the preliminary report [2], we shall, however, put 8 = 8 for a gas of constant 


specific heats. 











180 ERLING DAHLBERG {Vol. XIX, No. 3 


adsS 
S dz 


y — 
and the wall drag parameter 
g = B,a2f{/D > 0, 
where f and D are the conventional gas-dynamic expressions for the drag coefficient 
and the hydraulic diameter, F = pu’ S2f/D. (See e.g., A. H. Shapiro, Compressible fluid 


flow, New York, 1953, Art. 6.2) 
The differential equations for the normalized variables are found to be 


bh dD B, 
du* (us cal )(u* _ - . be + u*(gu* — ¥p*) 


a 


dx p* — (Bo — Bo/B)u* 
(2) 
b Bo b bs) ( Bo 
x + as e ” a * 
dy? (. »(p + 8 u b. B + ip" + 8 u* }(gu yp*) - 
” dx p* — (Bo — Bo/B)u* oe 
The expressions for the main physical quantities are 
RT = p*u*bo/Bo 
eae ; 
dh = 3 bo d(p*u*), 8 = B(p*u*). 
For 6 = const., h= p*u* bob Bo 
iES = ck’ S(u* — b/b,)b,/b = puS(u* — b/b,)bbo/aBo 
the electric energy extracted from unit mass per unit length being 
tES f{B d(p*u*) 1 du** 1 e 
pus 8, dz 2 dx dx 2 


M* = (8B, — Bo/B)u*/p*. 
For 8 = £B», we also have 
s = R In [(p*u*)’/p*| + R In S + const., 


ds aa (u* — b bo)” + gut 


a = *,,% ’ 
dx p*u 
‘ 8 
pu = Bp,/bop*(1 + u*/2p*)’, 
where p is the isentropic stagnation pressure. 


We note that the critical line (see Appendix) of the system of differential equations, 
Eqs. (2) and (3) 


p* = (Bo — Bo/B)u* 
corresponds to 


M*=1 or w = Bp/(B — 1)p 
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and that the jump relations, Eqs. (1), give 
Pil, = P22 , 
(pt — pi)/(us — ut) = —B, 

| (8/Bo) d(p*u*) = — (uk? — ux’), 

v1 
where subscripts 1 and 2 denote the state on the upstream and the downstream side 
of the discontinuity respectively. 
For 8 = const., the last relation may be replaced by 

(p* + p’) y (Bo sia Bo B)(u% + u*) 2 = 0, 

i.e. the jump is centered on the critical line. 


Because of the second law of thermodynamics, jumps may occur only from supersonic 
to subsonic motion 


Pi -~ is = Bo B,)ury < 0 < p> aoa (Bo aa Bo, B.)us ° 


For 8 = Bo, the change in entropy proves to be 
x Dati 
(s. — s,) = 2RB gee Tai (y"" — 1), 
p> 2n + 1 
where y = 8/(B 1) is the ratio of specific heats, and 


A = (pt — p*)/(pt + pir 
is a measure of the shock strength. 
We also note that 
(dp*/dx)/(du*/dx) — —B) when M—1. 
A search for singular points, characterized by 
dp*/dx = du*/dx = 0 


for the physically interesting cases (FE, B and o finite or zero) leads to the following 
results. Let 


G = p* — (By — Bo/B)u*. 
Then 
G(dp*/dx) = G(du*/dx) = 0 
gives either 


a) (0 < (b/bo)*/a < @) 


u* = b/d, 
vie 95 
a” ~ab 


and possibly some points on the line G = 0, which we shall call pseudo-singular points, 
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or 
b) ((b/bo)*/a = 0) 
view ot], 


—p* = 
a! a 


r. 

Phase diagrams. Numerical integration of Eqs. (2) and (3) would now yield the 
path of the phase point in the p*u*-plane for any particular flow situation with given 
initial values, ete. Without further simplification, little general knowledge is gained and 
the advantage of the normalized treatment is rather dubious. Because of the explicit 
x-dependence of the right-hand members of Eqs. (2) and (3), the phase point trajectory 
through any given point is not unique. However, in some cases it may be justifiable to 
neglect this z-dependence and we can then construct a phase diagram with unique 
trajectories. 

If, for instance, the ratio between width and height of the channel does not vary, i.e. 
Ay/Az = const., we can assume 


ES'” = const. 
BS'” = const. 

and consequently 
E/B = const. 


B*/pu = const. 


Thus 
@ = Goo/e 
b = bo 
(except in case of degeneracy, where we may have b = Oorb = o. If b = bh < 0, the 


physically relevant part of the phase diagram is the third quadrant, p* < 0, u* < 0. 
We shall, however, not discuss this case.) 

If, further, the area variation is moderate it may be permissible to use mean values 
for ¢g/a and y/a, especially as ¢ and y will mostly have the character of perturbation 
parameters. Finally, when computing ¢ and y, we must either assume that the electric 
conductivity is constant or approximate it by some function of p* and u* only. 

Let us assume that we may take the conductivity to be constant in this connection 
and also that the fluid may be treated as having constant specific heats (8 = 8) = const.). 

First, we consider the degenerate case, where the conductivity or the applied electric 
field is zero, which means that the power yield is identically zero. Equations (2) and (3) 
are then reduced to 


du* u*(u*¥(o + 1)/a — p*y/a) - 
eS = (2 } 

dx oy = 15 = 1)u* 

dp* (p* + u*)(u*(o + 1)/a — p*y/a) ‘es 
ee ee . ~ . . (0°) 

dx p* — (8 — l1)u* 


Here, , is an arbitrary constant and the relevant parameters, (g + 1)/a and y/a, are 
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finite. The equations can be integrated directly, using Eq. (4), which is reduced to 
p*tu* + u**/2 = const. (4’) 


(The integration constant may be put equal to unity because of the arbitrariness of by). 
Equation (4’) gives the shape of the trajectories, which are all similar. Evidently, if 
v/a = (¢ + 1)/a = O, every point is an equilibrium point, i.e. all flow variables are 
independent of x. Other possibilities are illustrated in Fig. 2, where the direction of 
motion along the trajectories is indicated by arrows. All diagrams in this and the follow- 
ing figures have been constructed for 8 = 8) = 4, corresponding to y = 1.33. 

A second degenerate case is characterized by finite conductivity and electric field 
but zero magnetic field. Equations (2) and (3) are now 


du*¥ C+ u*(u*g/a — p*y/a) 


ics : | Me bh | 9" 
dx p* — (6 — 1)u* ; oe") 
dp* _ _ BC + (p* + u*)(u*e/a — p*/a) (3””) 
dx p* — (8 — 1)u* . " 
(where C = cE’/pub; may be put equal to unity because of the arbitrariness of bp). 


The possible types of phase diagram are given in Fig. 3. 

Thirdly, we have collected in Fig. 4 the types of phase diagram that obtain when 
either or both of » and yW are so small as to be negligible. When both ¢ and y are zero 
(Fig. 4.2), an exact integral of Eqs. (2) and (3) can be found 


p*u* + u*’/2 — p*/B — u* = K = const. (5) 
or, solving for p*, 


_K+u* —«**/2 
pe u* — 1/8 


Substituting this in Eq. (2) with ¢ = y = 0 and integrating, we find 
[ ao dx 
J oa 


° 1 1 > 92 
oj E —4 , 6-34 Kp" t.. ne =e) | a“. 

/ u* — | (8 — 1) u* — 1/8 u (u* — 1/B)° 
It is evident from Eqs. (2) and (3) that, for small values of ¢ and y, the phase diagram 
will be substantially unaffected except near u* = b/b) = 1 (and, for y ¥ 0, for large 
p*-values combined with relatively small values of u*). Thus, except in special cases, 
it may often prove convenient to use Eqs. (5) and (6) for a first approximation*. On 
the other hand, it is obvious that the degeneracy, manifested by the line of singular 
points in Fig. 4.2, is removed by the introduction of a finite value of ¢ or y, no matter 
how small. 

Finally, when both ¢ and y are finite, the phase diagrams will be similar to those in 

Fig. 5. In this figure we have also included a sketch indicating what regions in 


(6) 


2In the earlier report [2] six diagrams are presented, which, for 8 = 4 and different values of K, 
give (B*/pu) fo dz, BRT /b? and M as functions of (h + u?/2)/b*. The diagrams also give (pu)ib/Bpsi 
as a function of (h; + u,?/2)/b*, where subscript 7 refers to the state at the duct inlet and p, is the 
isentropic stagnation pressure. 
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Fic. 2. Phase diagrams for degenerate cases, 1; co E =O. A phase diagram shows the variation along the 
channel of the normalized velocity, u*«< wu, and the normalized pressure, p* « pS, (where S is the 
cross-section area) for given values of the friction and channel divergence parameters, ¢ and y. In this 
figure we can recognize the normal gas-dynamic types of adiabatic flow. Thus, a represents frictionless 
convergent flow and ¢ constant area flow friction; in both cases the gas velocity tends to approach 
the speed of sound, represented by the critical line. Further, f represents divergent frictionless flow, 
accelerated when supersonic and decelerated when subsonic. In 6 the effects of friction and area change 
(convergent) cooperate, whereas in d and e they oppose each other. The singular line in d and e cor- 
responds to a balance between friction and area-change effects such that u* and p* are independent of z. 
0), the effect of the magnetic field, represented by 1/a in the 


With short-circuited electrodes (E 


parameter (y + 1)/a, is analogous to that of friction. 
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Fig. 3. Phase diagrams for degenerate cases, II; B = O. These diagrams illustrate the effect of Joule 
heating, o E?, in the absence of magnetic forces. In a this is the only effect and we see that the velocity 
approaches the speed of sound. The effects of other factors— friction, c convergent channel and d friction 
and convergent channel—which separately show the same tendency, will only modify the trajectories 





slightly. In e the effects of Joule heating and channel divergence oppose each other and we see that the 

former dominates at low pressures and velocities and the latter at high pressures and, less markedly, 

at high velocities. Friction, when small, will mainly influence the behaviour at low pressures and high 

velocities, f; its influence will gradually spread to higher pressures as the ratio between the friction and 
channel divergence parameters increases, g, h and 7. 
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Fia. 4. Phase diagrams; ¢ y = O. (N.B. Contracted p*-scale in d.) These diagrams represent the cases 
where friction and/or channel divergence are negligible. If these effects are small, the purely electro- 
magnetic effects illustrated in a will dominate except near u = E/B (u* = 1). This may be seen from 
b (friction), d (divergent channel) and f (convergent channel). Near u = E/B and when some other 
effect is great (c, e, g and h), the phase diagram will differ radically from a. One should note the smooth 
transition from supersonic to subsonic motion, which may occur with small friction, b, or small con- 
vergence (not shown). 
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Fic. 5. Phase diagrams; ¢ ¥ *O.(N.B. Contracted p*-scale.) These diagrams represent the cases where 

both friction and area change as well as electromagnetic effects must be taken into account. The values 

of » and y for the different diagrams in this and the preceding figure are indicated by the letters in the 

map of parametric regions. In t, k (r and s) we note the existence of a stable equilibrium point foru = E/B 

(u* and p* independent of z), in I, v (and ¢) it is unstable; this point corresponds to the singular line in 

Fig. 2, d and e. We also note the smooth transition from supersonic to subsonic motion that is charac- 
teristic for the regions I, t and a small part of region m (not illustrated). 
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the gy-plane correspond to the various types of phase diagram illustrated; the letters 
refer to the phase diagrams in this and the preceding figure. (The regions r, s, t, w and x 
have not been illustrated; they differ from their respective lower neighbours 7, k, 1, 
and q mainly in that the slope of the trajectories is negative along the whole u*-axis.) 
We note that the singular point (u* = 1, p* = ¢/y) is a stable nodal point in the super- 
sonic range with no limit cycles for parameter values in regions 7, k, r and s. When the 
parameters fall in the /, ¢ or v regions, it is a saddle point in the subsonic range and, 
for y < 0, it falls outside the physically relevant part of the phase plane. Regarding 
the pseudo-singular points on the sonic line, we find that no such points exist for the 
parametric regions n, g, v, w and x; when they exist they are two in number, and the 
one nearest the origin is always in the first quadrant (0 < u* < 2/(8 + 1)) and is always 
a “saddle point’’. The one farthest from the origin may be a ‘‘saddle point” (in regions 
k and s, u* > 1, and in regions 7 and r, u* < 0), a “nodal! point, stable from below,” 
(in regions / and ¢ and a small part of m near the ¢g-axis), a ‘‘focal point’’ (in the main 
part of region m) or a ‘‘nodal point, stable from above,”’ (in a very small portion of region 
m near the bottom left hand corner); for all the latter cases we find 2/(8 + 1) < u* <1. 

The pseudo-singular ‘‘nodal point’’ has a special significance in that the phase point 
moves across the critical line with a finite velocity, provided the initial point lies within 
a certain range. This implies the possibility of a continuous transition from supersonic 
to subsonic motion (or vice versa, depending on the type of ‘‘nodal point’’) without a 
throat in the duct and even when heat exchange effects are negligible. The parametric 
region where the pseudo-singular point is ‘“‘nodal and stable from above”’ is, however, 
so narrow as virtually to exclude the possibility of a transition from subsonic to super- 
sonic motion. (For 8 = 8 = 4, the width of this region is, approximately, 


—(8 — 1)/48 < y < —(6 — 1)/48 +7 10° 


for ¢ = O and decreases for increasing ¢.) Further, the uniqueness of a trajectory dis- 
appears when it crosses the critical line at a pseudo-nodal point; thus the phase diagram 
does not tell what will happen downstream from a smooth transition, but it gives some 
indication as to what can happen. Physically, this may be determined by end conditions, 
but there might also exist a corresponding real indeterminacy. 

In the general case, the phase point cannot cross the critical line except by dis- 
continuous. jumps from the supersonic to the subsonic range. This means that there may 
arise phenomena, somewhat analogous to the choking and the standing shock fronts 
known from compressible fluid flow in constant-area ducts with friction. (See e.g. 
Shapiro, loc. cit.) The analogy is far from perfect, however, as the electromagnetic 
“friction” force is negative for velocities smaller than F/B and also gives rise to energy 
exchange with the surroundings. A qualitative discussion of some typical cases of power- 
yielding, constant-area, small-friction flow, corresponding to the right-hand part of 
the phase diagram of Fig. 4.b, was given in the earlier report [2]. Keeping the boundary 
conditions fixed, we discussed how the phase-diagram image of the flow varies when 
the length of the electrode section (or, rather, the product of length and conductivity) 
is increased. 

For a duct fed through a convergent nozzle and discharged into a constant-pressure 
receiver, we showed that the flow will gradually be choked and the exit velocity de- 
crease until the velocity at the duct inlet corresponds to u* = 1. When the length is 
increased further, the flow rate will remain almost constant and the exit velocity will 
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increase, aS more and more energy is taken up by the gas from the external electric 
circuit. The existence of choking effects has been questioned by Neuringer [3] on the 
basis of his numerical results. Their absence there is, however, due to his assumption 
of constant inlet Mach number. 

If the gas enters via a convergent-divergent nozzle and is let out into a low-pressure 
receiver, the (supersonic) exit velocity will first decrease while the inlet velocity remains 
constant. Provided that the trajectory lies above the range of smooth transition of the 
pseudo-nodal point, a shock front will develop at the outlet when the velocity there 
becomes sonic and, as the length is increased further, this shock front will move up- 
stream. As long as u* > 1 everywhere in the duct, the outlet velocity, while remaining 
sonic, will at first increase slightly because of friction effects, but, if the shock enters 
the nozzle, the (still sonic) outlet velocity will again decrease. The shock may reach the 
nozzle throat and vanish there; then the flow will be choked in analogy with the con- 
vergent nozzle case. Depending on the area ratio of the nozzle and the initial state of 
the gas, it may also happen, however, that the shock stops in the divergent part of the 
nozzle or in the electrode section. This implies that the velocity in the duct inlet or on 
the downstream side of the shock corresponds to u* ~ 1 and, after this state has been 
reached, further increases in the length will not affect the flow upstream from this posi- 
tion, whereas the (sonic) outlet velocity will increase, as more and more energy is taken 
up by the gas from the external electric circuit. 

If, however, the initial trajectory lies below the separatrix of the pseudo-nodal 
point—which, for small values of g, is approximated by the no-friction trajectory 
through (u* |, p* = B — 1) with K = (8 — 2)(6 — 4)/8, Eq. (5)—no shock front 
should develop when the exit velocity becomes sonic*. Instead, the phase point should 
pass through the pseudo-nodal point with a finite velocity (a dp*/dr ~ ¢/(8 — 1)). 
Beyond the pseudo-nodal point it will continue along that trajectory which makes the 
outlet velocity sonic. Thus, when the length is increased further, the phase point will 
follow successively higher-lying trajectories and the (sonic) outlet velocity increase, 
as more and more energy is taken up by the gas from the external electric circuit. That a 
smooth transition is the only physically plausible phenomenon in this case is evident, 
as a jump from a point on the (supersonic) trajectory should end to the left of the vertical 
line through the pseudo-nodal point, where the trajectories go to infinity, and it would 
be impossible to satisfy the end condition of sonic outlet velocity. 

If the gas cannot be treated as having constant specific heats and/or constant con- 
ductivity, it is still feasible to construct a phase diagram for any given situation— 
provided that it is permissible to regard 6 and o as functions of p* and u* only. In most 
cases it should be sufficient to take into account the temperature dependence (7' « p*u*), 
but should it prove essential to include also a pressure dependence, it is necessary to 
disregard the area variation in this context by using p « p*/S > p*/S, . We have not 
included any such phase diagrams, as they would appear as distorted images of those 
presented, lacking in generality and, qualitatively, without any significant new features. 

Acknowledgments. I wish to thank Prof. H. Alfvén for his kind interest in the 
work, Dr. S. Lundquist and Prof. G. Borg for many helpful discussions and Mrs. K. 
Forsberg for drawing the figures. 

As the significance of the pseudo-singular points was not analysed in the preliminary report [2], 
this case was not included in the original survey and an erroneous statement was made (p. 14) to the 
effect that all trajectories in the region u* > 1 end on the critical line. 
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APPENDIX 
Remarks on a class of systems of non-linear differential equations. The system of 
non-linear differential equations studied in this report can be generalized as 


x’ 


II 


rie, y) R(a, ”)| (s ) 
1 


y’ = Q(z, y)/R(z, y)) 


(See e.g., Minorsky, Non-linear mechanics, Ann Arbor, 1947. Part IV. Chap. XX & 
XXII.) In addition to this system, we have the requirement—physically founded—that 
two functions, say F(x, y) and G(z, y), be continuous everywhere, whereas discontinuities 
in x and y may occur under certain, yet unspecified, circumstances. 
The phase trajectories of S, are the same as those of 
t! = Plt, D) i 
(Sz) 
y’ = Ua, y) 


but the direction and magnitude of the phase point velocity may be different for the 

two systems. The topological methods for sketching the trajectories of S, are well-known 

and straightforward. (See e.g., Minorsky, loc. cit. Part I.) Let us assume that we have 
made an outline of the phase diagram of S, and that it includes the isoclines P = 0 and 

Q = 0, and, further, that we have found and sketched the isolated points and curves 

where R = 0. By studying the sign variation of P, Q and R, we can now find the direction 

of motion of the phase point along the trajectories in the regions bounded by the curves 

P = 0,Q2 = Oand R = 0. 

Before examining the behaviour of the solutions at non-normal points, we shall 
make a number of assumptions, which shall hold in the whole of (the relevant part of) 
the zy-plane. 

A.l F, G, P, Q and R are continuous. 

A.2 The number of points where P = Q = 0 is finite. P and Q have continuous first 
derivatives at such points and P,Q, — P,Q, does not vanish there. 

A.3 R = 0 corresponds to a finite number of isolated points and a finite number of 
curves with a finite number of intersections. R possesses continuous first derivatives 
near points where R = O and continuous second derivatives, which do not all 
vanish, near points where R = R, = R, = 0. 

A.4 P,Q, R, R, and R, do not vanish simultaneously. 


. G. W. Sutton, The quasi-one dimensional flow of an electrically conducting gas for the generation of 
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A.5 At points where R = 0, the functions F’ = (F,P + F,Q)/R and G@’ = (G,P + G,Q)/R 
exist and are continuous. This implies sufficient differentiability of the functions 
involved. 

A.6 R, P, F, and G, do not vanish simultaneously. 

A.7 R, Q, F, and G, do not vanish simultaneously. 

We can now distinguish between four types of non-normal points of S, . 
I. A singular point is a point where P = Q = 0 but R ¥ 0. Its properties are similar 

to those of the corresponding point of S, , but its stability will depend on the sign of R. 

II. A critical point is a point where R = 0 but neither both R, and R, nor both 

P and Q vanish. A line of such points is called a critical line. The phase point velocity 
along a trajectory through a critical point changes sign (and is infinite) at that point; 
thus, the phase point cannot proceed continuously through (or stop at) a critical point. 
The assumption commonly employed is that the phase point jumps discontinuously 
from the critical point to some other point determined by the continuity of F and G. 
In our case, the jump occurs before the phase point reaches the critical line; this type 
of jump is physically plausible when distance and not time is the independent variable 
and the end conditions can influence the state immediately behind the jump (e.g., sub- 
sonic motion after the discontinuity). An interesting property of critical (and pseudo- 
critical, see below) points is evident from the assumption (A.5). Denoting a unit vector 
parallel to the trajectory through the point by s, we find 


s-grad F = s-gradG = 0 


which means that either /(G@) has a stationary value or grad F (grad G) is perpendicular 
to the trajectory. Should either of the gradients, say grad G, vanish at more than a 
finite number of points where R = 0, we may replace G by a suitable combination, such 
as H = F + G. Thus, normally the curves F = const. and G = const. (H = const.) 
are tangent to each other and to the trajectories at their intersection with curves corre- 
sponding to R = 0. If F or G is linear, say F = y — kx, R = 0 is obviously an isocline 
(dy/dx = k). We also see that jumps originating at critical points are impossible if no 
pair of curves, F = const.;G = const., possesses more than two points of intersection. 

III. A pseudo-singular (pseudo-nodal, etc.) point is a point where P = Q = R= 0 
but FR, and R, do not both vanish. Given the assumptions above, we can prove that 
any intersection between R = 0 and P = 0 or between R = 0 and Q = 0 is a pseudo- singu- 
lar point. If we assume e.g., R = P = 0, the assumption (A.5) gives F,Q = G,Q = 0; and 
Q must vanish, as, according to (A.6), F, and G, do not both vanish. The phase point 
velocity at the point will depend on the slope of the trajectory; using (A.2) we find that 
it will differ from zero, and using (A.4) that it will be finite, except if the trajectory 
should happen to be parallel to the critical line. In Fig. 6 we have illustrated the differ- 
ent types of pseudo-singular points, corresponding to the three types of singular points. 

a) Pseudo-focal (pseudo-vortex) point. The point is surrounded by arec-shaped 
trajectories beginning and ending on the critical line. 

b) Pseudo-nodal point. This point—like an ordinary nodal point—can be character- 
ized by two slopes: the “‘slope of approach” and the slope of the separatrix. We may 
expect a smooth transition if the phase point initially lies within the cross-hatched 
region. In a boundary-value problem like ours, the end condition may determine which 
trajectory the phase point will follow after passing the pseudo-nodal point; in an initial- 








192 ERLING DAHLBERG [Vol. XIX, No. 3 




















a) Pseudo-focal b) Pseudo-nodal c) Pseudo-saddle 
point point point 


Legend: ——-~——— Phase point trajectory 
ieee — Critical line 
altace: Line of characteristic slope 


Fig. 6. Pseudo-singular points. 


value problem this may be undetermined or determined by factors neglected in the 
analysis. 

c) Pseudo-saddle point. This point is again characterized by two slopes, correspond- 
ing to the two separatrices separating the regions where the trajectories end on the 
critical line, turn back, and begin on the critical line respectively. 

It is not to be expected that the separatrices, being discrete trajectories, correspond 
to physically observable smooth transitions. The application of restrictions having 
the nature of desiderata, such as constant gas velocity throughout the duct in our case, 
may easily hide the nature of the pseudo-singular point and lead to the presumably 
erroneous conclusion that a smooth transition will occur at a pseudo-saddle point (see 
e.g. Sutton [6]). If it should happen that the slope of the critical line equals one of the 
characteristic slopes, and especially the ‘‘slope of approach” of a pseudo-nodal point, 
further considerations may be necessary. One should, however, first of all ascertain 
that it does not depend on a never fully realizable desideratum, such as perfect. sym- 
metry of a circuit. 

IV. A pseudo-critical point is a point where R = R, = R, = O but neither all of 
R,, , R,, and R,, nor both P and Q vanish. It was proved above that R = PQ = 0 im- 
plies R = P = Q = 0; thus, neither P nor Q vanishes at a pseudo-critical point. Isolated 
points where R = 0 are pseudo-critical, as are points of intersection of critical lines; 
there may also exist lines of pseudo-critical points. Except in cases where the trajectory 
is parallel to a critical or pseudo-critical line, the phase point velocity along a trajectory 
through a pseudo-critical point, though becoming infinite at the point, does not change 
sign. The motion of the phase point is integrable and the pseudo-critical point will not 
affect. the continuity of the solution. In contrast with the discrete trajectories leading 
through a pseudo-saddle point, the discrete continuous trajectory leading through the 
pseudo-critical point at the intersection of two critical lines may well have a physically 
observable counterpart. The motion of the phase point along a trajectory near that 
through the pseudo-critical point should, in general, be quasi-continuous, i.e. when it 
reaches the first critical line the phase point will jump to a point on the same (S,)-tra- 
jectory beyond the second critical line (Fig. 7a). To prove this we assume that F% + 0 
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Fia. 7. Intersection of two critical lines. 


and that the trajectories are not parallel to either critical line. Denoting values at the 
pseudo-critical point by subscript 0, distance from that point by p and differentiation 
along the trajectory by d/ds, we find, using (A.5) 


dG/ds = (Gb/F¢) dF/ds + 6, 


where 6 is of the third order in p, whereas in general the difference between grad G and 
a constant times grad F is at least of the first order in p (second order in p if Fo = Fyo = 
Gio = Gy = 0). A typical plot of the variation in F and G along a trajectory traversing 
the two critical lines near the pseudo-critical point is given in Fig. 7b. (Geometrically, 
dF /ds approximates a hyperbolic paraboloid—such as dF/ds « xy—with its saddle 
point at the pseudo-critical point and its horizontal generatrices along the critical lines.) 
Thus, to a high order of approximation, jumps will occur between points on the same 
trajectory. In the first approximation, the length of the jump is 1.5 times the distance 
between the critical lines measured along the trajectory. Obviously, we cannot say that 
this jump must be unique, but in a physical application it presumably would be. 

When studying non-linear systems, it is essential to keep the physical background 
in mind. When a system is assumed to be stationary in time, the assumption may have 
to be tested experimentally. The justification of the discontinuity treatment of critical 
points is empirical and the same must be true regarding pseudo-singular and pseudo- 
critical points. As an example where pseudo-singular points occur in an initial-value 
problem we may choose the Abraham-Bloch multivibrator. (See e.g., Minorsky, loc. 
cit. Sec. 137.) Using the two grid voltages as phase variables, we can describe the situation 
thus. For vibrations to be possible, the characteristic parameter fia, = TRSmax/(R + 7) 
must exceed unity. When fa. is slightly greater than one we expect a pseudo-nodal 
point in the first quadrant with a smooth transition towards the origin and when finax 
is somewhat greater we expect a pseudo-nodal point with its slope of approach approxi- 
mately parallel to the critical line. Finally, when f,,, is still greater, this point corresponds 
to a pseudo-saddle point with a pseudo-focal point on each side. 
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Introduction to linear programming. By Walter W. Garvin. McGraw-Hill Book Co., 
Inc., New York, Toronto, London, 1960. xiv + 281 pp. $8.75. 


The book grew out of a course for selected personnel of the Standard Oil Company of California. 
It does not make great demands on the mathematical background of the reader and explains linear 
programming in elementary algebraic terms. 

Part I (82 pp.) deals with the general linear programming problem. The simplex method is developed 
and sensitivity analysis is discussed. An illustrative gasoline-blending problem is treated in considerable 
detail. Part II (65 pp.) is devoted to problems with the structure of the transportation problem, in 
which the triangularity of the basis leads to simplifications in the computational procedure. A generalized 
transportation problem with semi-triangular basis is formulated and a method for solving it is presented. 
Part III (118 pp.) treats a variety of special questions such as the presence of upper bounds, statistical 


and parametric linear programming, the revised simplex method, the resolution of degeneracy, and 





duality. 
W. PRAGER 


An introduction to the calculus of finite differences and difference equations. By Ixenneth 
S. Miller. Henry Holt & Co., New York, 1960. viii + 167 pp. $4.50. 


This book is not for one whose only interest is a working knowledge of finite differences for appli- 
cation to numerical computation. Rather, it is an elegant and rigorous presentation of the theory, 
developed with emphasis on the beautiful analogy with differential and integral calculus and differential 
equations. 

The calculus of finite differences is introduced in the first chapter, the theory of infinite products 
and its application to that calculus are demonstrated in the second, and Bernoulli polynomials and 
the Euler-Maclaurin formula treated in the third. The last chapter is devoted to linear difference 
equations in the real domain. 

Each chapter concludes with a large number of problems. Although written with all the pains a 
pure mathematician might demand, the book will appeal to applied mathematicians and numerical 
analysts who wish to gain a deeper knowledge of finite differences. 

WALTER FREIBERGER 


Bessel functions—zeros and associated values. Part III. Royal Society Mathematical 
Tables, Volume 7. Edited by F. W. J. Oliver. Cambridge University Press, New 
York, 1960. 1x + 77 pp. $9.50. 


It is well-known that the numerical determination of the zeros ot a rational polynomial presents 
difficulties of much higher order than the evaluation of the function at a set of values of the argument. 
This 1s even more true of transcendental functions like Bessel functions, and this volume therefore 
differs from its predecessors in a significant way. The 60-page introduction describes in fascinating 
detail how the difficulties inherent in this problem have been brilliantly overcome, and is well worth 
reading for its own sake as a stimulating essay in a baffling field of numerical analysis. 

For the practitioner, the tables will be invaluable. The main tables give the zeros of J, and of Yn 
with corresponding values of J/ and Y/; the zeros of J’ and Y/ with corresponding values of J, and Y,; 
and the zeros of the spherical Bessel functions~/ {x/2z} Jn41/2and+/ {w/2z} Yas1/2 with the corresponding 
values of their derivatives. 

The lay-out and typography are as exemplary as one has come to expect from the Cambridge 
University Press. 

WALTER FREIBERGER 
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ON THE THEORY OF PLANE STRESS* 


By 
EDWARD L. REISS AND STANLEY LOCKE** 


Institute of Mathematical Sciences, New York University, and Republic Aviation 


1. Introduction. The classical theory of plane stress [1] is concerned with a thin 
elastic plate subjected to edge forces which deform it so that there is no normal dis- 
placement of its middle plane. This theory is conventionally obtained from the exact 
three dimensional linear theory of elasticity for homogeneous and isotropic materials,’ 
by neglecting some of the compatibility equations and assuming that the shear and 
normal stresses, transverse to the middle plane of the plate, vanish and the remaining 
stresses are independent of z, the coordinate normal to this mid-plane. These assump- 
tions imply that only two boundary conditions, which are functions of one variable, 
are required along the edge. This is in contrast to the three of the exact theory which 
are functions of two variables. Therefore, a solution of the plane stress theory is not, 
in general, a solution of the exact theory.” The plane stress solutions can provide accu- 
rate approximations if the plate thickness is “‘small’’ compared to each of its lateral 
dimensions. 

In this paper, we indicate the relationship between the two theories’ by systemati- 
cally deriving from the exact theory the differential equations and boundary conditions 
of plane stress. Our technique of derivation also provides a systematic method of ob- 
taining increasingly accurate approximations to the exact stress distribution. This 
technique was previously developed and applied to other related elasticity problems 
(3, 4, 5]. In a sense, it is a generalization of the boundary layer method used by 
l'riedrichs [6] and later by Friedrichs and Dressler [7] in a study of the bending of plates. 
In the latter reference, it is also shown that by decomposing each stress component 
and edge force into the sum of an even and an odd function of z, the exact theory for 
the general problem of the loaded plate separates into two distinct systems of differ- 
ential equations and boundary conditions. One system describes the “bending” of 
plates and is partially treated in [3], [6] and [7]. The second system, which is concerned 
with the “extension”’ of plates, is considered in this paper. Hereafter, each stress com- 
ponent and edge force is assumed to possess its specific even or odd property. These 
are listed at the end of Sec. 2. 

Briefly, our method is to expand each stress component in a power series in the 
plate thickness, h. Substituting these expansions into the exact theory and equating 
coefficients of like powers of h yields a sequence of systems of differential equations to 
determine the expansion coefficients. The equations of the lowest order approximation 
(the ‘‘zeroth order interior problem"’) coincide with those of the plane stress theory. 
Higher order systems provide ‘“‘corrections’” to the plane stress approximation. In 
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Hereafter referred to as the exact theory. 
There are special cases for which the plane stress solution is a solution of the exact theory, e.g. [2]. 


In this connection see [9] and the discussion in [1]. 
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general, the expansions cannot satisfy the boundary conditions of the exact theory 
along the edge, nor do they represent the stresses in a region adjacent to the edge.* 
To determine, from the exact theory, appropriate boundary conditions for the previously 
mentioned differential equations and to obtain approximations for the stresses near 
the edge, we employ a second expansion in h, whose coefficients are now functions of 
new independent variables. These new variables are obtained by “stretching” the co- 
ordinate normal to the edge of the plate.*° Substituting these expansions into the exact 
theory and equating coefficients of like powers of h results in a sequence of boundary 
value problems, whose solutions yield approximations to the stress distribution near 
the edge. The requirement of existence of single valued solutions of these problems 
provides the desired boundary conditions. The boundary conditions obtained from the 
lowest order approximation coincide with those of the plane stress theory. 

In Sec. 6, the procedure for obtaining the higher order approximations for the three- 
dimensional stress distribution is summarized. 

2. Formulation. An elastic plate of thickness 2h is considered as a three-dimensional 


elastic body bounded by the planes z = + A and the cylindrical surface x X(s), 
y = Y(s). The bounding planes z = handz = —h are referred to as the faces of the 
plate. The surface, | z| < h,x = X(s), y = Y(s), is called the edge, and the ‘‘smooth” 
curve B:z = 0,2 X(s), y = Y(s), is called the boundary curve. 


For simplicity, it is assumed that the faces are force free’ and that the edge is sub- 
jected to an arbitrary equilibrium distribution of normal and shear forces. We further 
assume that the magnitudes of the edge forces are sufficiently small to insure that the 
exact theory remains valid. 

If we introduce the dimensionless variable 


oh? M) 
then the faces are given by ¢ = +1. If the stresses are functions of x, y and ¢ the equa- 
tions of the exact theory are [1] 

Equilibrium equations, 
G.2.. + hl{c,.. + ¢.,.,| = 0, O.y.¢ thle... + o,.,] = 0, (2) 
O,¢ t+ hlo.ss + Oy1) = 0; | 
Compatibility equations, 
O22, t h[Ao, + 2,,,] = 0, o,.+ +t h[Ac, + 2,,,] = 0, 
O2,.er + Qiee + h? Ac, = O, Cee.t¢ t+ AQ». + h? Ac,, = 0, (3) 


Ory is¢ + hQ ty + h’ Ao,, as 0, Try,tt ai h*[Ac,, + | wi 0. 
Here we have employed the notation, 


a” oO 
A= 73 t ; Q = 


y4 ° 2 he, + o, + ¢,),; 
Ox oY 


Vv 


+\/- 


4See Secs. 3 and 4. 
5 See Sec. 4. 


°This restriction is not essential, but merely condenses the calculations. 





1961] ON THE THEORY OF PLANE STRESS 197 


and an independent variable appearing as a subscript and following a comma indicates 
partial differentiation with respect to that variable, e.g., ¢,.,.; = 00,,/0¢. 

To complete the formulation we require appropriate boundary conditions. These 
are obtained by specifying the applied forces as 


g(x,y, + 1;h) = o,(2z, y, &1;h) = o,,(z, y, & 1;h) = 0; (4a) 
o,(X, Y,¢;h) = LOe, ¢;h), o,(\, Y,¢;h) = Qs, ¢;h), 
aja, 7,8 R(s, €; h). 


(4b) 


In (4b), » signifies the outward normal to B and s the distance along B. L, Q and R 
are the prescribed normal and shear edge loads. Since Eqs. (2-4) are the stress formulation 
of the elastic problem for the plate, no expticit reference to Hooke’s law is required. 


“ce 


In this paper attention is focused on the “extension”’ of the plate by the edge forces. 
This implies [7] that L, Q, o, , o, , 0, , and o,, are even functions ef ¢, while R, ¢,, and 
, are odd functions of ¢. In the ensuing analysis we shall make frequent use of the even 
or odd property of each of these quantities without explicit reference. 
3. The interior problems. We assume that each stress component symbolized by, 


o(x, y, &; h), can be represented, in an asymptotic sense, by a power series in h 
a(x, y, (> h) ~ > a(x, y, Oh. (5) 
Furthermore we define o' 0 if i < 0. The functions o'(x, y, ¢) are called the interior 


stress coefficients and are assumed to possess the same even or odd property as 
s(x, y, €; h). We also assume that the prescribed edge forces in (4) can be expanded in 
power series In A: 


Ls, c;h)| (L'(s, a) 
126 5h) = > Q'(s, ph’. (6) 
UR(s, ¢; h); R'(s, $)) 


Substituting Eqs. (5) into the differential equations (2) and (3) and boundary con- 
ditions (4a) and equating coefficients of the same powers of h yields, from the coeffici- 
ents of h’, 


? n—} a—1 n n—-1 n-1 
0:2 T Oz 2 + Fas ee 0, Orzy,t + Ory,z i Cyy = 0, (7a b c) 
As , J 


> 


yett +. Ao? 7 + ieee = 0, 


-| 
be 
Q 
+ 
~ 
a 
an 
Qq 


(8a, b, ¢) 


, , , n—2 r n- n-2 
te 2 Q) a + Ao., = 0, Ory tt oa o + Ao,, = 0, 


c tr + Ac: + 9g -— 0; 


a(t, ¥, £1) = oo, (a, y, 1) = oz, y, +1) = O. (9) 


(8d, e, f) 


We now proceed to systematically analyze (7-9). Since o° = 0 if i < 0, it follows 
from (7) and (9) that, 


Cie = Oy = 6, = 6, = 0. (10) 
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An immediate consequence of (8a, b, f) and (10) is that, 


o; T Six, y), oy i Si(x, y); Ory = Bust; y), (1 1) 
(1 + y)Q'(x, y) = Siz, y) + Sita, y), if n=0,1. 
Equations (7a, b), (9) and (11) then imply that, 
1 ~ 1 a ) 
w.. = 6, = @. (12) 
We can, in fact, prove’ that, 
4 ; ‘ l »» 2 9 
C.2 = Oy =—6,= 9), Q = Oz, y) = ——(g. + g,), n=0,1,°--. (13) 
l+yp 


Hence (7c) and (8c, d, e) are identically satisfied by (13) for all n and (7a, b) become, 


n 


o> Cas = Gas Tr Co = 9, n=0Q0,I1,°---. (14a) 


Zy,2z ' wv 


Furthermore if o% and o” satisfy (8a, b) then using (13), we see that, 


AY =0, 2n=0,1,-:-. (14b) 
Integration of (8a, b, f) with respect to ¢ yields 
: e 

- af et’ 
a” : a? ee S.(2, 7 = 'o Q’ ig _ | | Ao" “es UF B) dB at. 

Pe oF pt’ 
o(z,y, f) = S,@,y) - o 0 | | Ao) (x, y, 8) dB dg’, (15) 
n i ee 
Got. ee Dealt 7/) = ma OF as. ) m= QQ, 1, «++, 


and (13) as solutions of (7-9) provided that S” , S* and S?, satisfy the following differ- 
ential equations obtained from (14): 


n 4 c 8 me a a 2 ee 
See te — (O°, — oO"), | + | | Ao. a6 at’ . 
2 Ox P x 3 ; 
ef ee ee , (16) 
Pe) = ae ss Bay eas [Q — vi .| + Ao’; dp de’, 
2 Oy ' o 4 : ; 
AM" = A(S°+ S) =0, n=0,1, 
In particular, we observe that (15) reduces to (11) for n = 0, 1 and (16) reduces to 
S2. + S2,., = Sty.2 + Si, = AVS? + S})) =0, if n=0,1. (17) 


Equations (10), (11) and (17) give the equations of the classical theory of plane stress 


[1]. It is convenient to refer to (17) with n = 0 and n = 1 as the differential equations 
of the zeroth and first order interior problems respectively. By recursive application 
of (15) and (16) we obtain all higher order interior problems. For example, if n = 2, 3 
o, = Sz, y) — 1(AS?~* + allo sg 
o, = Siz, y) — 1(AS*~* + Q' os ? 
, 
Co S?Ar y) + ; 2 sil 


iSee App ndix. 








1961] ON THE THEORY OF PLANE STRESS 199 


S; z + ae = Dosis + By. _ A(S? + S)) = 0. 


4. Formulation of the boundary layer problem. In the previous section we obtained, 
for all n, explicit differential equations, (16), and expressions for the interior stress 
coefficients (13) and (15). These results were derived from the exact theory independent 
of the boundary conditions along the edge, (4b). For each n, the differential equations 
form a “fourth” order system and hence require two boundary conditions on B, which 
are functions of only s. The exact theory requires three conditions, (4b), that are func- 
tions of s and ¢. Furthermore we observe that the interior stress coefficients have a 
specific dependence on ¢, (13) and (15), while the prescribed stresses in (4b) are arbitrary 
functions. Therefore the interior expansions (5) cannot, in general, satisfy the edge 
boundary conditions. In fact, if the series (5) represent the solution they do so in some 
region away from the edge, i.e. in the ‘interior’ of the plate. Near the edge, i.e. in the 
‘boundary layer,”’ the solutions vary rapidly from those of the interior solution to 
satisfy the boundary conditions of the exact theory. To obtain an expansion that uni- 
formly represents the solution up to and including the edge we employ a “stretching” 
of the independent variables. Our procedure, which has previously been employed in 
studying other related problems [3, 4, 5], differs somewhat from that of Friedrichs [6, 7]. 
We assume that the rapid variations that occur in the boundary layer depend only on 
the direction normal to the edge. Therefore, only the independent variable in this di- 
rection is stretched.” 

To facilitate this, consider any point P on B and assume, without loss of generality, 
that the x and y axes coincide, respectively, with the normal and tangent to B at P. 
We define at P the boundary layer coordinates ~, y, ¢ by “‘stretching’’ the variable x 
so that, 


br 
=> (18) 


Then for sufficiently small h, any fixed neighborhood of the edge, no matter how small, 
in the x variable corresponds to an arbitrarily large one in the £ variable. Roughly this 
implies that the boundary layer effects penetrate from the edge into the plate a “dis- 
tance”’ of order of magnitude h. 


If f(g, y, &; h) = o(a, y, £; h) is the generic symbol for the stresses as functions of 
£, y, ¢, we assume that 
fEy~65h)~ Lie y, on. (19) 


Here f' are called the boundary layer stress coefficients and we define f*' = 0 if ¢ < 0. 

Introducing the transformation (18) and the expansions (6) and (19) into the exact 
theory (2-4) and equating coefficients of like powers of h, we find that the coefficients 
of h” yield 


7 f + frie ™ —f; “ ’ rs + ee si —fy ’ ) + Fiat _ —fy } (20a, b, c) 


Vitle= 2, Vh=-f2—-Ts, Vi+ln=—f2 (lab, 

co72y7 ! yn pn—2 25n yn-1 n—2 

V i r= —Jis.sv 3 “on “< ~ Fae ‘ 

Ves I fg J v3 Vi I ig ’ (21d, e, f) 
ge = a vine tn ; 


8For more detailed descriptions of these methods see [3-8]. 
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filé, yy 1) = fig, y, +1) = fig, y, & 1) = 0; (22a) 
f:(0,0, ) = L'(0, 9), f:.(0,0, §) = R’O, 5), f:,(0,0, 9) = QO, 9), (22b) 
where, 


: a” f a°f - l n n n 
Vf = 5+ Pr and, M99 = Ts,5 i. +h > f)- 


OF 
Equations (20-22) are written for y = 0, i.e. in the normal plane to B at P. 

To complete the formulation of the boundary layer problem we require that each 
f approaches or ‘‘matches” the corresponding interior stress as h — 0, or equivalently 
from (18) as § — —o. More precisely, the matching conditions, or the asymptotic 
forms of the boundary layer stress coefficients are obtained [3, 4, 5] by assuming that 
each interior stress coefficient, o'(x, 0, £), possesses a Taylor series expansion about x = 0 


o(z,0,8) = >, Aida’, (23) 
where, 
fs. oft 
1 0’o'(0, 0, &) , , 
7 = , >0 | 
A;(¢t) = + Ox’ » Hite (24a) 
0, : $<. 
It follows from (5), (18), (23) and (24a) that in the neighborhood of x = 0, 
o(z,0, ¢;h) ~ > DE, 0, Or’, (25) 
where 
Dé, 0,9 = >> AT (HE (24b) 
are the interior stress coefficients near x = y = 0 as functions of (é, ¢). If we define the 
“reduced boundary layer stress coefficients” F"(£, 0, £) as, 
F°(é,0,6 =f'é,0, 0 — D’é, 0, 5, m2 @, 1, +++, (26a) 


then employing (19), (25) and the definitions of f(£, 0, ¢; h) and o(x, 0, ¢; h) the asymp- 
totic form for the f” is obtained by identifying each f" with the corresponding D” near 
x = 0ash—O. We write this condition, using (18) and (26a), as: 
lim F’(é, 0, ¢) = 0. (26b) 
ose 
In obtaining (26) we assume that the terms in f” which vanish as & ~ — © do so faster 
than any negative power of &. 
The boundary layer problem, (20-22), (24) and (26), is systematically analyzed 
in the next section. 
5. Analysis of the boundary layer. For each n, the boundary layer problem separates 
into two distinct systems which are called Problem P" and Problem 7”. Problem P" 
involves only the coefficients f* , f” , f" , and f”, and consists of Eqs. (20a, b), (21a, b, ce, d), 
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the first two of (22a, b) and (24) and (26) for these coefficients. Problem 7” is concerned 
with the two remaining coefficients and consists of (20c), (21e, f), the last of (22a, b) 
and (24) and (26) for these coefficients. 

We shall associate with the F” of Problem P” a function ¢"(é, 0, ¢) which may be 
the solution of the following boundary value problem on the semi-infinite strip D, 
Is|<1,é<0 

V'‘¢" _ 0; 
@(€,0, +1) = hk’), ¢ er(€,0, + 1) = 0; lim [6"s: , d"e] = 0; (27) 
f--@ 


o";:(0, 0,¢) = a’(0, $), ¢":(0, 0, g) ” B'(0, g), 


where a” and 8” are respectively even and odd functions of ¢ and {> k"(€) dé = 0. It can 
be shown by an elementary calculation that if ¢",, , ¢";, and ¢", are uniformly continuous 
functions of ¢ and if ¢", ¢", and ¢", are single valued functions, then, 


1 
[ a"(0, ¢) dt = 0. (28) 
v-—1 
Using (11), (13) and (17) we can show that 
Fr =@1; Fl = @e » Fi. = —@ite » F) = Ve" (29a) 
is a solution of Problem P° if ¢° is a solution of (27) with, 
k(é) = 0, a’(0, ¢) = L'(0, §) — S20, 0), B°(0, ¢) = —R*(0, 0). (29b) 


Application of (28) to (29b) gives the first boundary condition at y = 0 for the zeroth 
order interior problem as, 


S20, 0) = ; [ L°(O, §) dé. (30a) 
Similarly, by employing (11), (13) and (17) Problem 7° has the solution, 
Fy = V's Foy = We (31a) 
if y’ is a solution of the following boundary value problem on D 
Vy’ = 0; 
vix(E,0, 1) = 0; lim {oe , Wir} = 0; (31b) 


I] 


v°.(0,0, 9 = Q°O, § — S20, 0). 


Existence of a solution of (31b) requires that the integral around the boundary of the 
normal derivative of ¥° vanishes. This yields the second and final boundary condition 
for the zeroth order interior problem at y = 0 as, 


vO , 1 ; 0 
$°.0,0) = 5 | Q%0, ak. (30b) 


The solution to Problem P' can be reduced to solving (27) with 
kg) = p°(é, 0, + 1), a'(0, % = L'(0, ) — Si0, 0) + v°,0, 0, o, (32a) 
B'(0, s) = —R'O, 8) 
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and 
__ 1 oe 0 nb = 1 a 0 AR! =e e 
F, — rr Vi» ] F, Q zz Voy ] Fu Pr ] (32b) 
F, = »V'¢' + 2y’, . 


From (28) and the second of (32a) we find the first boundary condition at y = 0 for 
the first order interior problem as, 


1 al 
S:0,0) =5] (’'O,s) + 0,0, Hl de. (33a) 
a J-1 
If we set 
Pec = zit esi vo. ty ? pe bg Oe ae vty ’ (34a) 


- ° ° 0 ml ee ml se af fa . . ° 
then (34a) is a solution of Problem 7" if Z'(é, 0, ¢) is a solution of the following boundary 


value problem on D: 


ll 


V7E! = (1 —») l a 


='.(é,0, +1) = 0, lim (ae, ar} = 0, (34b) 
='.0,0, 6 = —Q'(0, ) + S!,(0, 0) +» | R°O, ¢) dé. 


ov =§ 


It can be shown by an elementary calculation that if =', is a uniformly continuous 
function of ¢ and if Z is a single valued function then, 
al 
vl / ] l 1 ai 79 « ' _ 99 
S.,(0,0) = 5 (Q'(0, §) + rfR0, $)] de. (33b) 


) 


“4 J—i 


This is the second and final boundary condition for the first order interior problem. 
Higher order boundary layer approximations may be obtained by proceeding in a 
similar manner. 

From Problems P” and 7” given above we observe that, in general, F, , F,, and F,, 
do not vanish, while in the Appendix it is shown that the corresponding interior stresses, 
g, , o-, and a,, do vanish. Thus there is a contribution to the transverse normal and 
shear stresses only in the boundary layer and we would therefore expect the plane 
stress theory, i.e. the zeroth order interior problem, to yield a ‘‘good’’ approximation 
to the exact solution in the interior. This has been verified for the problem of an infinite 
plate with a circular hole stretched by a uniaxial force at infinity. The three-dimensional 


corrections to the plane stress theory yield only “‘small’’ changes in the stress concen- 
tration factor. These results will be reported in detail elsewhere. 

6. Summary. In the preceding sections we have shown, using systematic expansion 
procedures, the relation of the plane stress theory to the exact theory. The plane stress 
theory appears as the zeroth order interior problem, (11), (13) and (17) with n = 0 
and (30), the solutions of which supply a first approximation to the three-dimensional 
stress distribution. If we introduce an Airy stress function, G’(z, y), such that, 


¥0 yO 0 0 0 0 
S, = Ces ? S, = 7 ’ Bos = —G a ? 


then the zeroth order interior problem can be written in the familiar form 
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A’G® = 0; 
and 


10 yO ] ‘ 0 % vO 0 ] : 0 
S, = pee = 9 [ L df, Szy = —G..y = di Q df, on B. 


In addition, the expansion procedure provides “‘three-dimensional corrections” to 
the plane stress theory. The first correction is in the boundary layer and is obtained 
from the solutions of P° and 7°. The next correction is obtained from the solutions of 
the first order interior problem, (11), (13) and (17) with n = 1 and (33), and Problems 
P' and T". We notice from (32) and (34) that the solutions of the zeroth order problems 
provide inhomogeneous terms for the first order problems. We define an Nth approxi- 
mation to the exact theory by introducing o“’(x, y, ¢; h) as the generic symbol for the 
stress components of this approximation where, 

N 


o(r,y, 5h) = DY [o(a, y, ) + F(a/h, y, Oh". 


n=0 
Here o"(x, y, €) and F"(x/h, y, ¢) are, respectively, the interior and reduced boundary 
layer stress coefficients of order n obtained from the solutions of the nth order interior 
and boundary layer problems respectively. 
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APPENDIX 


Recalling that o” and Q" are even functions of ¢ and that o°%, and o7, are odd func- 


tions of ¢, we prove the Theorem. o°, = 02, = 0) = Oand Q" = Q"(a, y) for all n. 

The proof is by induction. Let o7, = o7, of = os’ = o:*' = 0, and M = Q(z, y) 
and Q"*’ Q”*'(x, y). Then from (7c) and (9) we immediately conclude that 
( il 0 and hence from (8c), 2"** = 2"**(x, y). It follows from (8d, e) and (9) 
and the inductive assumption that, 03? = o%3* = 0. The induction is completed by 


reference to qs. (10-12). Thus, 


rae > oh = 0, a lead yp > o;,h* itd 0, od died pS o:h* = 0. 


i=0 1-0 +=0 
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BOOK REVIEWS 


(Continued from p. 194) 


The control of multivariable systems. By Mihajlo D. Mesarovic. Technology Press of 
M. I. T., Cambridge, Mass., and John Wiley & Sons, Inc., New York and London, 
1960. xi + 112 pp. $3.50. 


The monograph is divided into two parts. The first part contains some general philosophical dis- 
cussion culminating in the observation that control of a complex system can be materially improved 
if more is known of the internal structure. To illustrate this point, the author considers several different 
types of structures together with their transfer functions and control policies. 

The second part is devoted to the synthesis problem. This is first formulated as the problem of 
minimizing a functional of a collection of state variables, and then is simplified to that of minimizing a 
quadratic functional. A number of specific calculations are then carried out. 

There is no doubt that the author has some interesting ideas in this modern fascinating field of 
control theory. However, the book is so loosely organized and the results so imprecisely stated that it 
is exceedingly difficult to discern what they are. This is a book which required a strong editorial hand, 
and did not receive it. 

RicHarD BELLMAN 


Wave propagation and group velocity. By Leon Brillouin. Academic Press, New York 
and London, 1960. xi + 154 pp. $6.00. 


The purpose of this book is to bring to the fore some of the fundamental old papers on the subject 
of wave propagation. The title is perhaps misleadingly short, since the book deals with the distinction 
and relations between phase velocity, group velocity and signal velocity of a propagating wave. The 
author presents the subject against the background of its historical developments, as he is indeed well 
qualified to do. 

The structure of the book is somewhat unusual. The first chapter gives several basic definitions 
concerning wave propagation, followed by a summary of the theoretical situation around 1910. Particu- 
lar emphasis is placed on the difficulties that were encountered, in those days, in making a clear dis- 
tinction between group velocity and signal velocity in dispersive media. The following four chapters 
cover the fundamental work in this area by Sommerfeld and by Brillouin, which culminated in a clari- 
fication of the above concepts. This part is presented in the form of four papers, originally published 
between 1914 and 1932; these papers are reprinted with very little editorial comment. Finally, the last 
chapter discusses briefly the application of the ideas developed in these papers to guided waves (both 
acoustic and electromagnetic). 

The work on this subject by Summerfeld and by Brillouin is, of course, a classic in the field. 
Presenting it in this form, under one cover, will undoubtedly facilitate the task of those who are interested 
in the complete, original treatment. 

C. ELBaum 


Weight-strength analysis of aircraft structures. By F. R. Shanley. Dover Publications, 
Inc., New York, 1960, xiii + 404 pp. $2.45. 


This is an unabridged reprint of the first edition (McGraw-Hill, 1952), to which a new preface and 
bibliographies on optimum design of structures and creep buckling have been added. 


(Continued on p. 220) 
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A THEORY FOR THE SQUARE GUARDED HOT PLATE—A SOLUTION OF THE 
HEAT CONDUCTION EQUATION FOR A TWO LAYER SYSTEM* 


By 
I. G. DONALDSON 
Dominion Physical Laboratory, Lower Hutt, New Zealand 


Summary. A theoretical treatment is given for the action of the guarded hot plate 
apparatus used for measuring conductivities. The flow of heat through the metal hot 
plate itself, and the effects of temperature unbalance between hot plate and guard plate 
are both taken into account. The temperature at any point in the sample or in the heater 
plate may be determined. The theory is comprehensive and any of the dimensions or 
controlling parameters may be varied in the calculation. 

Introduction. Although the thermal conductivity of insulating materials has been 
determined for many years with the guarded hot plate apparatus, it was not until recently 
that any attempts have been made to assess theoretically the errors associated with 
the use of this apparatus. Without the benefit of theoretical assessments the standards 
that have been specified by the ASTM [1], RILEM [2] and the British Standards Insti- 
tution [3] have from necessity been drawn up relatively arbitrarily. Thus they differ to 
some extent one from another and in the light of recent theoretical and experimental 
work, do not appear to be sufficiently stringent for the 1% accuracy often demanded 
of the guarded hot plate apparatus. 

Figure 1 illustrates the layout of a typical “guarded hot plate 
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Fic. 1. Schematic diagram of the guarded hot plate apparatus based on the design of the apparatus at 
Dominion Physical Laboratory. The z scale has been exaggerated (XX 3) to show the construction more 
clearly. 
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slabs of the material under test are placed one on either side of an electrically heated 
metal plate. Outside again are two metal plates kept cold by a flow of water. Thermo- 
couples measure the temperature difference between the hot and cold plates, and the 
conductivity of the material is calculated from the measured temperature difference, 
from the known rate of generation of heat per unit area of the hot plate and from the 
thickness of the specimen. 

In this calculation the thermal gradient and the flow of heat are presumed to be 
everywhere parallel to the axis of the instrument (perpendicular to the plates) but in 
fact some lateral flow of heat occurs near the outer edges of the hot plate and of the 
specimen slabs. Consequently the hot plate is provided with a guard ring in imitation 
of electrical practice. This is a plate in the form of a hollow square closely surrounding 
the hot plate, fitted with an electric heater and kept at the same temperature as the 
hot plate. In effect the hot plate becomes merely the central portion of a larger heated 
plate and the transverse heat flow in its vicinity is much reduced. Of course, an attempt 
is made to prevent heat passing through the edge faces of the specimens by surrounding 
the whole apparatus with thermal ‘insulation’, but the thermal conductivity of such 
material is appreciable. Indeed thermal measurements tend to be much less accurate 
than the equivalent electrical measurements, mainly because the range of thermal 
conductivities of available materials is much less than the range of electrical conduc- 
tivities. 

Earlier experimental tests. An experimental check by Robinson and Watson [4] 
involving 20 different guarded hot plates, all built to the ASTM specifications, found 
that 75% of the measured conductivities were in agreement to within 3% of the mean 
but some of the results were out by as much as +13% and —16%. Gilbo [5] showed 
that small temperature differences between hot plate and the guard plate (which will 
be called ‘‘unbalance’’) could cause considerable errors; for his guarded hot plate 0.2°F 
temperature difference caused 3% error. He also found that the ambient temperature 
around the apparatus was another large source of error. Roux et al. [6] investigated 
errors due to temperature unbalance of the order of 0.03°F and reported that the error 
increased as the conductivity of the sample decreased; for cork (thermal conductivity 
0.30 B.t.u. in/hr ft?°F) this small unbalance caused an error of more than 0.5°%. Pascal 
[7] in tests on cork 10 em thick in his 50 em X 50 cm hot plate with gap filled 
with araldite found that an unbalance equal to 1% of the temperature difference be- 
tween hot and cold plates caused an error in conductivity of 35%. Woodside and Wilson 
[8] investigated unbalance effects and studied the dependence of these errors upon the 
magnitude and direction of the unbalance, the size and design of the heater plate, the 
conductivity and thickness of the specimens tested; and the temperature difference 
between hot and cold plates. Errors as high as 6% were found under conditions that 
satisfy the ASTM requirements. 


Earlier theoretical assessments. 
1. Edge heat loss. 


Theoretical assessments that have been made may be divided into two groups; 
those predicting errors due to edge heat loss (due that is to heat flow through the edge 
insulation) and those determining the effects of temperature unbalance. Somers and 
Cyphers [9] did a full analytical treatment of the problem of edge heat losses for a par- 
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ticular case, where the temperature of the edge of the sample was equal to that of the 
cold plate. Dusinberre [10] used a relaxation technique for a similar problem, but allowed 
for some edge insulation. Pascal [7] also used relaxation to determine the minimum 
value of the ratio of guard ring width to specimen thickness which would make the 
transverse flow of heat, negligible in the test section of the specimens, but his treatment 
is two-dimensional not three-dimensional so there is some doubt about the validity 
of the results. 

In the treatment by Woodside [11] the cold plate is taken to be at zero temperature, 
the hot plate to be at a uniform temperature 6, , and the edge faces of the sample are 
taken to be at a uniform intermediate temperature e@, . These assumptions enable a 
two-dimensional solution to be found by a Schwarz-Christoffel transformation. The 
solution for the particular case when e is zero agrees well with that of Somers 
and Cyphers [9], but the assumption of a uniform temperature on the edge face of the 
specimen seems unrealistic. 


2. Unbalance effects. 


The effect of a temperature difference between the hot plate and its guard ring 
has been computed by Pascal [7] from the lateral flow of heat through the gap between 
them. This, however, neglected the lateral flow through the specimens themselves. 

Woodside [12] on the other hand based his theory for the error due to unbalance on 
the simplified system shown in Fig. 2 which takes account of the specimen but ignores 
the flow of heat across the gap. He justified the extension of the specimens and plate to 
infinity on each side of the gap and the assumption of two-dimensional heat flow by 
relaxation calculations (not given in his paper). By applying two successive Schwarz- 
Christoffel transformations, Fig. 2 was transformed to a simpler system for which the 
solution was readily obtained. As an experimental check the change in the total heat 
flow due to a temperature difference between hot plate and guard plate was measured 
on three sets of apparatus and was found to agree within 5° with the predictions of 
the foregoing theory. 

With the exception of the results of Pascal which are in doubt for the reasons given 

bove, in all the above theoretical approaches it has been assumed that the temperature 
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Fic. 2. The artificial plate-specimen temperature system assumed by Woodside [12] to determine un- 
balance effects. The solution obtained for a uniform temperature over the entire guard plate—hot plate 
combination (@ = 6, say) with cold plate maintained at @ = 0 is added to this solution. The reduced 


boundary conditions and dimensions are shown. 
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of the hot plate and of the guard plate are uniform over the face of each plate. This 
assumption is unjustifiable, especially in cases where the guard plate temperature is 
different from the hot plate temperature or where the temperature maintained outside 
the edge insulation of the apparatus differs from that of the hot plate, even if some 
edge insulation is used. To overcome these doubts it was felt that some account should 
be taken of the heat flow in the metal plate itself. It was found that in so doing it was 
possible to develop a theory which took into account both edge and unbalance effects 
and also matched very closely the method by which most guarded hot plates are 
controlled. 

Theoretical analysis. The general arrangement of a guarded hot plate apparatus 
has already been described in Fig. 1. During operation the cold plates are maintained 
at some fixed temperature and the heat generated electrically in the hot plate and guard 
plate is adjusted so that the temperature difference between the central areas of the 
hot and cold faces of the sample is maintained constant, while the temperature difference 
between corresponding points on the outside edge of the hot plate and on the inside 
edge of the guard plate is kept as nearly zero as possible. A theory has been designed 
to fit this practical model, with but one modification. It has been necessary to assume 
that the air gap in the metal plates is actually filled with metal, i.e., the plate heater 
and guard heater are sandwiched between two continuous plates. The effect of this 
will be discussed later. 

Because of the symmetry of the apparatus it is only necessary to consider half of it. 
The theory therefore discusses the double slab illustrated in Fig. 3, consisting of a metal 
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Fic. 3. A cross sectional view of the block used as the basis for the theory. ABDC is the sample and 
CDFE is the metal plate. 


plate of thickness ¢ and conductivity K,, and the sample of thickness ¢ and conductivity 
K, . Each has a length and width 2a. The origin of coordinate is taken to be the centre 
point of the interface. 

It is assumed that the free face of the specimen is at a constant temperature 7, 
while heat enters the face of the metal plate at a prescribed rate Q(x, y) which depends 
on position. Around the sides of the compound slab there is assumed to be an insulating 
layer of conductivity K, and thickness p whose outside face is kept at some temperature 
taken to be zero for the theory. All other temperatures are measured in degrees Fahren 
heit from this value. 

In both the sample and the metal plate the temperature obeys Laplace’s equation 
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VT = 0. (1) 


Boundary conditions have been prescribed for the upper and lower faces. At the 
interface the temperature and the normal component of heat flow must be continuous. 
An approximation will be introduced in considering the edge faces. It will be assumed 
that the flow of heat across these boundaries is proportional to the difference in temper- 
ature between the edge and the temperature outside the insulant. This approximation 
is discussed in a later section. The parameter controlling this will be termed the co- 
efficient of edge heat transfer and, writing h for this, the boundary conditions are 


T=T, z=, —a<z2z <a, —a<y<a, (2) 
aT 
K a = —Q(z, y) z= —, —-a<2x<a, —-a<y<a, (3) 
Cz 
oT —_ 
5 + hT = 0 x = +a, —a<y<a, —§<z< f, (4) 
Ox 
oT ~ 
a ~ AT = 0 z=-a, -a<y<a, —-§<2<f, (5) 
On 
oT on" § 
5, + AT = 0 y =a, —a<xz<a, —-§<z2<f, (6) 
oT ‘ine 7) 
ay = 0 y = —a, —a<2<a, -—E<2z2<f, (7) 
ol 
Tasca: = Teameic 2=0, —a<x<a, -a<y<a,_ (8) 


= QO, —a<2 <a, —a<y< 4. (9) 


rn 
| 


> her . (oT 
(7) = x(*) 
Oz metal Oz sample 


Q(x, y) is the heat input at the metal face and is defined as follows 


—b<2z<b, —b<y<b, Q(x, y) = 
) 
te < a 
= y<e 
—a<y< —c| Oz, ¥) = 4, 


cx <n 


—c<z<-t 
‘ SS . 9 == < y < Cc, Q(x, Yy) = 0 
»<s <6 
—a<y< —c} Qr.y =H 
cxiy< ea 
=—@G <. 2 Ss, “4 a y ce. Q(x, ¥} = @ - (10) 
e<e<etea J 


For either the metal sheet or the sample, the solution of Eq. (1) under the boundary 
conditions 4 to 7 is of the form (see for instance, Carslaw and Jaeger [13]) 
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T= > > (M,, cosh l,,z + N,, sinh l,,z) cos a,x cos B.y, (11) 
r=1 s=1 


where the coefficients 1 and N will be different in the sample and in the metal sheet 
and where 


V,=ar+ Be (12) 
and a, and 8, are the roots of 

atanaa =h, (13) 

6 tan Ba = h. (14) 


It is seen that a, = 6, . Different symbols have been chosen for clarity. In order to 
determine values of the coefficient (symbolised as M and N) fitting the boundary con- 
ditions noted in Eqs. (2), (3), (8) and (9) it is convenient to note (Carslaw and Jaeger 
[13]) that when a, is def heed in the above manner, the function f(x) obeys the trans- 
formation 








= 2(a2 + h’) cosa,t f° ,, . 
iz) = ; = l xX) COS a,x dx. (15 

’ 2, (a2 + h*)a +h J, 1 ) 
The boundary condition (2) states that at the surface z = ¢ the temperature 7 equals 
T, so that by applying Eq. (15), in both the z- and y-directions, this constant temper- 
ature can be written as the function 


=~ < ft (ae +h’) (8° + h’) 


OS a,% COs BY ‘ 
— sin a,a sin 6,a. (16) 


ra 
or ns eB {(a? + h?)a + h}f (8? + h? ja +h} 


Similarly at the surface z = —é the heat input is Q(z, y) defined in detail in Eq. (10) 
and this may in turn be written 





Ro = 





> — Ala? + h')(6; + h') cos d,x cos By _ 
ai — 0,8. {(a2 + h®)a + hh i (B; + h*)ja+h} ' 


+ q,(sin a,a sin B,a — sin a,c sin B,c)}. 


{qo sin a,b sin B,b 


Q(z, y) = (17) 


Now at any point within the sample, the temperature distribution can be written [see 
Eq. (11)] as 


o 


” oe 2 (A A,, cosh l,,z + B,, sinh l,,z) cos a,x cos B.y (18) 


=1 s=] 
but at the boundary z = ¢ the temperature follows Eq. (16), consequently 


= toler + h*)(6; + h*) sin a,a sin B,a__ , 
4,, cosh 1,,¢ + B,, sinh l,,¢ aBA(a? + ha + hi + Wa + hi} (19) 





Again, at any point within the metal plate 
= } p> (C,, cosh l,,z + D,, sinh l,,z) cos a,x cos B,y. (20) 
r=l1 #¢= 


The value of Q at the boundary z = —é can be deduced from this, using (3) and com- 
parison with expansion (17) shows that 
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—C,, sinh l,,£ + D,, cosh l,,é 





© _ Aa? + h’)(8) + h’) {qo sin a,b sin 6,b + gi(sin a,a sin B,a — sin a,c sin B.c)} (21) 
K,l,.0,8.{(a2 + ha + hb} + Wa + h} 
Further, the boundary conditions (8) and (9) applied to Eqs. (18) and (20) give 
Biss = Ges , (22) 


| pp. = K.B,, ° 


Relations (19), (21), and (22) allow the coefficients A, B, C and D to be expressed in 
terms of the athe constants, so that within the metal sheet 


p> 


‘a B.{(a, + h’)a + hy (f+ h ja + h} \(K,, sinh l.,é sinh 1,,¢ + K, cosh l,,é cosh l,,¢) 


e 


A(a? + h’ (8° + h’) cosa,x cos By 








] {qo sin a,b sin B,b + q,(sin a,a sin B,a — sin a,csin B,c)} sinh l,,¢ 
(23) 


+ K,T, sin a,asin B,a cosh Ls | cosh I,,z 
. ‘ , : : _ 
+]— KL. (ge sin a,b sin B,b + ¢q,(sin a,a sin B,a — sin a,csin B,c)} cosh l,,¢ 
+ K,T, sin a,a sin B,a sinh Lt | sinh Lah 
and within the sample 
P= 2. 2 
4(a, + h')(B, te) ) cos a,x cos B,Y 
-a,B,{(a? + h*)a + h}{(B? + h®)a + A}(K,, sinh l,,é sinh U,,¢ + ‘K, cosh l, .£ cosh l, 2) 


if ; ere 
( , {qo sin a,b sin B,b + q,(sin a,a sin 8,a — sin a,c sin B,c)} sinh l,,¢ 
\ ben 


(24) 
+ K,T,sina,asin B,a cosh | cosh I,,z 
. cere. ' ; . 
+ | —7~ {go sin a,b sin 8,b + q,(sin a,a sin 8,a — sin a,csin 8,c)} cosh l,,§ 
+ K,,7 sin a,a sin B,a sinh Ls | sinh La}. 
Thus for any point (x, y, 2) in the sample or in the metal plate, the temperature is given by 
T = u(x, y,2)q0 + v(x, y, 2) + Ma, y,2)To , (25) 
where u(x, y, 2), v(x, y, 2,) and X(z, y, 2) are, for any particular guarded hot plate func- 


tions of position only. In operating the apparatus, as pointed out earlier, the heat inputs 
into the hot plate and the guard plates are controlled to maintain a fixed temperature 
difference between the central areas of the hot and cold faces of the sample. For the 
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theory, as it is here, central points only are used. The theory may however be modified 
for any point or set of points. At the same time the temperature difference between two 
points, one on the edge of the guard plate and one on the edge of the hot plate, is main- 
tained zero. Applying these conditions to Eq. (25), two equations 


T, + AT = p20, 0, O)go + (0, 0, O)g, + ACO, 0, O)T> (26) 


and 


{u(c, d, 0) — u(b, d, O)}qo + {v(c, d, 0) — v(b, d, 0)}¢q, (97) 
af 


+ {r(c, d, 0) — A(b, d, 0)} = O 


are obtained and hence it is possible to determine g) and q, for any T, . Variation of 
T. is effectively the same as varying the temperature of the controlled atmosphere 
surrounding the apparatus. The measured thermal conductivity of the sample is then 
obtained from the heat input/sq. foot of the heater, the dimensions of the apparatus, 
the thickness of the sample and the temperature difference across the sample 
(K.,) = Goo 5 (28) 
= 1(b +c)” AT ” 

Unbalance effects. The theory as so far described covers all variations provided the 
apparatus is in balance, i.e., no temperature difference between a pair of points on the 
outside edge of the hot plate and the inside of the guard plate. In practice a gap fitted 
with a poor conductor (in most cases—air) separates the guard plate from the hot plate, 
but in the theory above no such gap exists. 

To match the theory more closely with the apparatus as it is in practice, it 
is necessary when determining unbalance effects to consider the problem in two parts: 
(a) the effect of the correct heat flux in the metal plate and (b) the correction for the 
effect of the temperature difference across the gap on the sample itself. 

(a) Thus, if 7’ is the temperature difference across the actual gap (at the pair of 
points) and the gap is filled with material of conductivity AK’, then the heat flows would 
be the same as if the temperature difference between the same two points in the solid 


metal plate were 


T’ =T"’ =. (29) 
Kn 
For the determination of the measured thermal conductivity in this case T’’ replaces 
the zero in the right hand side of Eq. (27) and this is combined with Eq. (26) to determine 
a qo . This g thus takes lateral heat flux in the metal plate into account but must be 
corrected for effects in the sample. 
(b) In the sample itself the actual temperature difference 7’ applies. Consider the 
sample alone with the boundary conditions (4)-(7) as before, but with 


7 = ~2=¢ —a<z<a —a<y<a, (30) 
T = T(z, y) z=0 —a<2x<a —a<y<a, (31) 


where T(x, y) is zero if | xz | and | y | are both less than 3(b + c) otherwise it is equal to 
T’. (T’ is positive when guard plate is the hotter.) 
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The solution obtained is 


Bn {(oe + h®)a + h}{(B2 + h®)a + A} sinh U,,.¢ (32) 
X {sin a,a sin B,,a — sin a,(b + c) sin 36,,(b + o)}. 


= 1T’(a. + h’)(B; + h’) cos a,x cos B,,y sinh l,,,.(¢ — 2) 


The total heat entering the sample from the hot plate must therefore be 


ov + K, [ [ j (223) dz dy (33) 
J—(b+e)/2 J-(b+e)/2 \ 02 J 2=0 


and the measured conductivity will be 


1¢ ‘is 
(K.)mess = ae ap '° % 
= 16K .T'lnn(@m +h )(B, +h’) sin $a,(b + c) sin 38,(b + o) woth 1..t (34) 
a a Be {lar + h’)a + hi {(Bz + h’)a + h} , a4 


< {sin a,a sin B,,a — sin }a,(b + c) sin 38,(b + c)}. 


The coefficient of edge heat transfer. In a system in which the insulation is in 
contact with only one other material a simplification often used is to assume that the 
insulation may be replaced by an extremely thin layer of material of almost zero con- 
ductivity, i.e., the flow of heat through the insulation is taken to be perpendicular to 
the interface between insulation and the material being insulated. In this case the 
boundary condition used pertains to this interface and is taken in the form 

dT 


dn 


= h(T — T,), (35) 


where 7 is in the direction normal to this interface and 7’, is the temperature at which 
the outside of the insulation is maintained. A reasonable approximation for the co- 
efficient of edge heat transfer is then 

hm ot (36) 

Kd 

provided that there is good contact between the two materials, and that the temperature 
gradient through the insulation is much greater perpendicular to the interface than 
that parallel to the interface. In Eq. (36) K, is the thermal conductivity of the insula- 
tion, K is the conductivity of the other material and d is the thickness of the insulating 
layer. 

In the guarded hot plate the insulating layer is, however, in contact with both the 
sample in the apparatus and with the metal of the guard plate. As well as this the temper- 
ature gradient parallel to the interface will be comparable with that through the insula- 
tion, thus any attempt at defining h as in Eq. (36) could lead to considerable inaccuracy. 
Nevertheless it is necessary to obtain a uniform value of h over the entire edge match- 
ing as closely as possible the edge heat losses for any particular set up in order that the 
equations may be solved. This value of h was determined as follows. The arrangement 
was simplified by working in two dimensions only. The sample/hot plate block LM M’‘L’ 
(Fig. 4) and the insulation @M’VV’ were considered separately. It was assumed that 
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Fic. 4. The boundary conditions assumed for the determination of the coefficient of edge heat transfer, h. 


the block LMM’L’ was edged by a thin layer with coefficient of edge heat transfer 
h in place of the layer of insulation and using this the heat flow over any section of the 
face MM’ was determined. The temperature at the point N was also determined. The 
block of insulation was then considered and it was assumed that the heat flowing over 
any section of the face 1/.M’ of this block matched that flowing out over the same section 
of the face MM’ of block LMM’'L’. The temperature at the point N’ was then determined 
and this temperature was equated to the temperature determined for the point N of 
the block LMM’L’. 

For the block LMRS (Fig. 4) the boundary conditions are (assuming that the co- 
efficient of edge heat transfer is h) 


at +=0 T=fT,, (37) 
, (aT fer 
at x= f¢ K,(2) = K(2 : (38) 
: VG? an Ox/ , 
ab z=¢ ( (39) 
+ 
at w= eC+E K, — = Q, (40) 
. Ox 
vi 
at y= -—a : — h?' = 0, (41) 
oY 
AT , 
at y=a - +hT = 0, (42) 
OY 


where K,, is the conductivity of the metal plate, K, is the conductivity of the sample 
and the subscripts m and s refer to the metal and the sample respectively. The heat 
input Q is so chosen that if the edge insulation were perfect the temperature difference 
across the sample would be 40°F. Similar boundary conditions apply for the 
block L’M’RS. 

Thus solving Laplace’s equation in the block LMM’L’ with these boundary con- 
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ditions the temperature at any point within the block is given by 


T= > {A, cosha(z — §) + X= B, sinh a(x — £)} cosa,y 0O<2*<¢ 


8 


= > {A, cosha,(x — §) + B, sinha,(x — §)} cosa,y c<r<seteé 


r=] 


y {A, cosha,x — ¢ 


r=1 


2t) — B, sinha,z — ¢+ 2t)} COS a, (43) 
c+t§EStrS St 2 


2¢) — a B, sinha,(x — t + 28)} COS a@,¥ 


+, 


> {A, cosh a,(x — 


r=1 


¢+2<2< 2¢ + 8), 


where 


an(K.7, cosh a,£ + o sah a.) 
A, = os a (44) 


‘a? + ha + h} cosa,a(K, cosh a,f¢ cosh a,f + K,, sinh a,f sinh a,é) 








on(—K.T. sinh a,§ oe ra cosh a.t) 
B, = wr , 2 ; —— = m A, (45) 


‘(a2 + ha + h} cosa,a(K, cosh a,f¢ cosh a,§ + K,, sinh a,f sinh a,£ 





and by differentiating this with respect to y and multiplying by the applicable value 
of the conductivity the heat flow over any section of the boundary MM’ may be obtained 


K@(%) = > K.a, sin auth A, cosh a,(x — §) + a B, sina(x — »} 


0 Uy 


0O<2r<¢ 


= > K,,a, sina,a{A, cosha,(x — §) + B, sinha,(x — §)} 
r=1 


rFszsctté 





(46) 
= > K,a, sin a,a{A, cosha(« — ¢ + 28) — B,sinha,(x — ¢ + 28)} 
ftiEszezssc tsa 
oY —— ; 
= y K,a, sin a,a) A, cosh a,(x — ¢ + 2) — K B, sinh a,(x — ¢ + 2) 
f+ 2%<2 < A+ 
while the temperature at N(y = a, x = £) is 
Ty = > A, cosa,a. (47) 
For the edge insulation MVV’M’ it is assumed that 
at y=a+p T=0 (48) 


at zxz=0 P - To (-y+a+p) (49) 
x= 2¢ + & . 
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while at 


a i 
y=a x,(*“) = K(a)(*") (50) 
OY OY 


i ‘ y=a ) 


where K(x) (07/dy),-. is defined above in Eq. (46) and the subscript 7 pertains to the 
insulation. 

The temperature assumed along the edges MV and M’V’ of this block, as given in 
(49) is from necessity an approximation and will give for any simple arrangement. of 
insulation a slightly low value of h. If the entire hot plate assembly is covered by the 
same thickness of insulation the temperature variation across these faces does not differ 
greatly from the linear variation assumed. A relaxation technique was used to determine 
this temperature variation in several cases. 

Solving Laplace’s equation in the block 1J/VV’M’ under the boundary conditions 
(48), (49) and (50) the temperature at any point within the block is found 


1 


T's ; 
T= (—y + a + p) 
p 
> > } 1(¢ + fa, sin a,a 
wae eet lnaK ; [n° + tas(¢ + £)°} 


" a ee ee , " 
[ax ma(¢ + §) sin > (A, sinh a,é + B, cosh af) + nr(1 — cos nz) 
= (51) 


5 ‘ ; . ‘ mre 
(KA, cosh a,f + K,,B, sinha,¢ — (K,, — K,)A, COs 57, uf )| 


“Al -- &)) 
onl {na(a +p—y)\. Nx 
. Ree sinn ‘ So Saari Sa a7. a 
, 46 + &T ol | 2@e+é) | 2(¢ + &) 
nr'p nrp 
I cosh I 
2(¢ —+ ¢) 
2 ré 


ry 


and thus the temperature at the point N’(y = a,x = £) is 


; ss ~4(¢ + )7T, nrp nre 
T=T,- > # —- = tanh ; = L a Oars a 
=— n3p 2A + &) oe T=) 
I >" .< 1(¢ + Ela sin @,a 
Ke: ( (nal[ne + 4ae(¢ + &) (BO) 
; . nr : , 
[ax “a(¢ + £) sin > (A, sinh af + B, cosa,£) + nr(1 — cos nm) 


. oe ; “ nro 
(-K A, ecosha,¢ + K,B, sinha,¢ — (K,, — K,)A 008 = ) 
(p | 


| -\ 
AS “TW eh 


nTp , nrc 
tal : 
X tanh + Mag +E 


As N and N’ are in fact the same point the temperatures as defined by (47) and (52) 
should be the same and thus by equating these h may be computed in terms of the 
other constants. This equation is however not linear in h, whereas it is linear in K,; hence 
it is much more convenient to choose an h and then solve for K; . The variation of K;, 
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with h was determined for various discrete values of the dimensions and for various 
materials as heater plate and sample. In all cases it was assumed that the apparatus 
was 1 foot square, i.e., a = 0.5 foot. 

The results of these computations show that provided the thickness of insulation 
is greater than the thickness of one sample (p > ¢), Ah does not vary with insulation 
thickness, neither does it vary significantly with the conductivity of the sample itself 
for the range of values considered (h < .04 ft™', 02 < K, < .10 Btu/ft°Fhr, .01 < K; < 
04 Btu/ft°Fhr K,, = 0(10°)K,). It was found that h is approximately proportional 
to the conductivity of the insulation K,; , inversely proportional to the conductivity 
of the metal plate K,, , and for the Q used here it is inversely proportional to T, + 40 
for the range —35 < T, < O°F. 

The variation of h with ¢ and ¢ is rather more complicated but is given by 

h = o(f, &), (53) 
where $(¢, £) is given for various discrete values of ¢, 1/16" < & < 3/8” in Fig. 5. This 
may all be combined in the single formula 


_ _Ko(s.8) 
h= KT, + 40) ”" 


which defines h to within 7% of the computed values. 
To justify to some degree the method of determining the value of h, K; has been 
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Fic. 5. The variation of ¢ (¢, &) with ¢ for discrete values of ¢. 
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determined by equating the temperatures at several other points at the interface be- 
tween the two blocks instead of at N and N’. The values obtained have been compared 
with those obtained above and it was found that for all the values of K; likely to be met 
in practice, i.e., 0.01 < K; < 0.04, Btu/ft°F/hr. the variation of h across the edge of the 
sample or of the metal plate for a particular K, is less than 6% of the value determined 
above. As in the simplified heating system used in this section (no guard heater) a vari- 
ation of h by 10% merely alters the temperature difference across the centre of the 
sample from hot to cold plate by less than 0.3%, it is to be expected that in the calcula- 
tions for the full apparatus a variation of h of the order of 10% will have a negligible 
effect and hence there is little justification for any attempt to determine h to a higher 
degree of accuracy. 

Conclusions. The theory developed above thus effectively matches most aspects of 
the guarded hot plate and its operation. It is possible to allow for the variation of any 
of the following 


(a) The outer dimensions of the guard plate. 

(b) The inner dimensions of the guard plate heater. 

(c) The outer dimensions of the hot plate heater. 

(d) The thickness of the metal plate. 

(e) The material of the metal plate. 

(f) The thickness of the sample. 

(g) The material of the sample. 

(h) The insulation around the apparatus. 

(i) The temperature difference between hot and cold plates. 

(j) The temperature of the cold plate (this effectively alters the temperature of 
the surroundings). 

(k) The positions of the thermocouples controlling balance. 

(1) The amount of unbalance. 

(m) The gap width and the material with which the gap is filled. 


This theory was not intended to be used in the direct measurement of thermal con- 
ductivity, but rather as a full check of the accuracy of the application of the simple 
linear flow formula under all conditions. The results will be used in the designing of an 
apparatus in which the total error will be minimised. 

The double series (24) is rather slowly coverging and preliminary computations 
showed that although 100 terms were sufficient to give the temperature to 4 decimal 
places, insufficient accuracy was obtained to permit gq and g, to be computed with 
any certainty. 

It is intended to obtain a full solution for a 1 foot square guarded hot plate on a 
digital computer, but at the present time access to one is being awaited. Once a full 
solution has been obtained for an apparatus of one size, the solution will be equally 
applicable to guarded hot plates of other sizes by scaling the results. 
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Stochastic processes: problems and solutions. By Lajos Takacs. Methuen & Co. Ltd., 
London, and John Wiley & Sons, Inc., New York, 1960. xi + 137 pp. $2.75 


A true description of this book is given in the Introduction as follows: “This book --- is primarily 
a collection of problems and their solutions, and is intended for readers who are already familiar with 
probability theory. Its aim is to summarise the fundamental notions and theorems of stochastic processes. 
The proofs of the theorems are generally omitted or only a brief outline is given. --- The scope of this 
book extends over the theory of Markov chains, Markov processes, stationary stochastic processes, 
recurrent processes and secondary stochastic processes. The problems are taken from the fields of 
natural sciences, engineering and industry.”’ 

This is an ideal book for anyone who is already familiar with elementary probability theory and 
would like a brief introduction to the above topics and their applications. By avoiding difficult topics 
and proofs, the author manages to scan an amazing range of subjects for such a small book. The treat- 
ment of Markov chains and processes, though highly condensed is very well done but the subject of 
stationary stochastic processes is perhaps too brief to be very instructive. 

G. F. NEWELL 


Tables of In T(z) for complex argument. By A. A. Abramov. Translated from the Russian 
by D. G. Fry. Pergamon Press, New York. 331 pp. $17.50 


These are six-figure tables of In I(x + ty) for x = 1 (.01) 2, y = 0(.01) 4, that were compiled at 
the Institute of Precision Mechanics and Computer Technology of the U. 8. S. R. Academy of Sciences. 
Formulas for computing /n I'(z) outside this rectangle are given. A loose sheet furnished with the book 


contains a nomogram facilitating interpolation. 


Technische Schwingungslehre. By Karl Klotter. Springer-Verlag, Berlin, Gottingen, 
Heidelberg, 1960. xv + 483 pp. $14.95. 


This second volume of a two-volume series on vibrations deals with vibrations of multiple-degree- 
of-freedom systems. It is intended as a textbook rather than as a handbook, so the arrangement of 
material is primarily pedagogic; moreover, emphasis is on fundamental ideas rather than applications, 
and certain more sophisticated topics such as Laplace transforms or non-linear differential equations 
are deliberately omitted. The topics are: Part I: meaning of degrees of freedom, with many examples; 
electrical-mechanical analogies; equations of motion, and their matrix form; Lagrange’s equations and 
Hamilton’s principle, with numerous examples; normal modes, with graphical methods—first for two 
degrees of freedom, and subsequently more generally; properties of linear differential equations; Routh 
stability criterion, and examples; quadratic forms, and small displacements; damping; practical examples. 
Part II: torsional vibrations and shafts; difference methods; iteration; graphical methods; bounds; 
matric algebra; more complex problems; literature survey. 

From the student’s viewpoint this is an exhaustive and exhausting book (almost 500 pages). It 
contains an interminable number of examples, many of which are useful and instructive (the book is 
particularly rich in graphical methods—with the surprising omission of phase plane methods, however), 
but which are excessive in number. For textbook use the reviewer would prefer a concise text on 
methodology with a reasonable number of selected examples; also, the omission of non-linear systems 
seems unwise. 

C. E. Pearson 


(Continued on p. 230) 
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AN IMPLICIT, NUMERICAL METHOD FOR SOLVING TWO-DIMENSIONAL 
TIME-DEPENDENT DIFFUSION PROBLEMS* 


BY 
THOMAS A. OLIPHANT 


Los Alamos Scientific Laboratory, University of California 
Los Alamos, New Mexico 


1. Introduction. We develop an implicit scheme for the numerical solution of 
two-dimensional, time-dependent diffusion problems and other related problems. The 
underlying principle is the solution of a general set of linear, simultaneous equations 
which may occur in nine-point differencing schemes. 

We apply the method to both linear and non-linear problems. Each time step is 
solved by an iterative procedure and a convergence condition is derived for the linear case. 

The stability of any time-dependent treatment is investigated in terms of its par- 
ticular differencing scheme independently of our method of solving the simultaneous 
equations. Our scheme is quite useful when applied to simple, unconditionally stable 
differencing schemes. Schemes involving more than nine points can be constructed and 
the method is in principle generalizable to three dimensions. 

2. Derivation of the method. The basic partial differential equation including the 
linear and non-linear cases which we wish to consider is 


00 
at 


V-{fr, V(@)} + Ar, 006- SH =A K (1) 
This is not the most general type of equation to which this method is applicable, but 
it illustrates the generality of the method. We are specifically interested here in ap- 
proximating the solution to this sort of equation by the use of finite differencing schemes. 
Generally speaking, implicit schemes are written in the form of a set of linear or non- 
linear simultaneous equations for the dependent variables over the mesh at an advanced 
time. Various implicit methods of this type have proved to be of considerable practical 
use. One such method is the Douglas-Peaceman alternating direction method [1]. It 
represents the application of the one-dimensional Bruce, Peaceman, Rachford, and Rice 
method [2] to two dimensions. In this method one direction is treated exactly at each 
time step. An alternative method constructed by Baker and Oliphant [3] permits an 
exact treatment over the entire two-dimensional mesh for a single time step. However, 
this method has a practical limitation in the symmetry which must be imposed on the 
second partial derivative term in order that the factorization which is involved can be 
carried out. The method which we propose here does not have this limitation. How- 
ever, as we shall see, this method is intimately related to the Baker-Oliphant method. 

We wish to approximate the partial differential equation (1) by a finite differencing 
scheme. We write our differencing scheme as a set of simultaneous equations 


Unless otherwise specified, we sum over repeated indices. The (4, 1) and (7, j) indices 
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refer to the position in two-dimensional space of the mesh points and the n index refers 
to the advanced time. The D,, contains such things as source terms and 67;* and 67;’, 
the values of @,; at earlier times. In non-linear problems, the coefficients Bj} and D,, 
may also depend on the 67; . In order to avoid confusion in the rest of this section, we 
will leave the advanced time index n understood. We will write out specific differencing 
schemes in Sec. 5. The present method is a straightforward generalization of the Bruce, 
Peaceman, Rachford, and Rice method to two dimensions. We do not factor the quantity 
Bii(Bii = Aj}B?) as was done by Baker and Oliphant. However, we must still use a 
nine-point differencing scheme rather than a five-point one in order to make the equations 
all consistent with each other. Therefore, in Eq. (2), 7 and j differ from k and I, re- 
spectively, by as much as, but no more than, unity. 

The first step in the solution of (2) is to triangulate the quantity Bj) , i.e., to write it as 


Bui = witbnn » (3) 
where w7? has all zero elements except when 
m=k,k—1 
(4) 
n=I1,l-—1 
and b,*? has all zero elements except when 
i=m,m+i1 > 
(5) 
j=n, ntl. 


Just as in the one-dimensional case, we set the diagonal element bj; (no sum) equal 
to unity (although other normalizations could be used). Next, we write D,, as 


Dar = Wiridmn « (6) 
Substituting (3) and (6) into (2), we get 
Wi DmnI:; = WEYmn » (7) 
so that we must have 
Gnn = DmnBi; « (8) 


Equation (3) can now be written out in algebraic form as a set of equations with 
no sums over indices 





wi = By; (9) 
wig) = Bit — wip bij | (10) 
wis? = By — wie bai 5 (11) 

wii = Bi — wey Bsa — wis bh — wi bY | (12) 
#+1j+1 

pis — Bah — way bi , (14) 

if ‘7 ’ } 


Wi; 
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Biit —_ wi, ipsit! 
pit) = Ae (15) 
W; j 
By = wae bi; (16) 
By = wi be (17) 


Equations (9) to (12) serve to determine w{;*’~', wij"’, wii7’ and wi! . Equations (13) 
to (15) serve to determine bif!’*’, bi}™, and bii*?. Unfortunately, however, Eqs. (16) 
and (17) involve only quantities which have already been determined by calculations 
at previous mesh points, and in general they are not satisfied as identities. This in- 
consistency can be removed if we go to an iterative scheme. To do this we write out 
Eq. (2) and add and subtract quantities from the left side. The resulting equation (18) 
is so far exactly equivalent to the original equation (2). 
—  - + Bi; lig ie 4 e — ry 1j+1 + e oi ee + Bi 1g Rew 

+ Biji0:; + Biy O41 + (Ba — Ba (18) 

+ tata ee + Bi Seats + ) aaa 941541 = D;; ° 
Transposing appropriate terms to the right hand side, we get 
Bip O.-aj-1 HF Bip O15 HF BiG Osta Big "Oss + BiiOis + Bir Ose 

+ Bi; i O41; -1 + Bij Osos + Bir Osarser — Ei; ’ (19) 


where 
Ey; = Dy — (Big ** — Bie) Oia + (BY — BE) Oa}. (20) 
Now, Eqs. (9) to (15) are as before, but instead of (16) and (17) we have 
ijt om wis, (21) 
eit elf“. (22) 


These equations determine §';"’** and ;}"’ 

Then, we solve a single time step for a linear problem in the following way. The 
values of @ appearing on the right hand side of (19) are now taken to be guessed values 
to be used in the iterations. After computing the 6’s and guessing the @’s, E,; of (20) 
is determined. In Eq. (6), replacing D,; by E,; , we obtain 





, ij-l1 ,t~t7 e~tf—8 
E;; = W;; Jij-1 at Ww; Ji-13 _ W:; Ji-1j-1 23 
iii stn Fi ; (2 ) 
Wi; 


With Z;, known, we can compute g,; throughout the mesh. Then, with the b’s and g’s 
known throughout the mesh, we can compute the 6’s by Eq. (8) which is written out 
as follows: 


. ij+l1 i+1; i+1j+1 ‘ 
0; ii bres bj; 05541 — bi; i413 oat bi; Oi eiien ° (24) 


The values of @,,; determined in this way can be used to find new guessed values. We 
will discuss the convergence of the iterative scheme for linear problems in Sec. 4. 

It is important to notice here that the Bj; are restricted only in forming a nine- 
point scheme, but are otherwise completely general. Therefore, there is considerable 
flexibility in writing appropriate differencing schemes. 
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3. The boundary conditions. Very general boundary conditions can be set up. 
We can set up our formulas to take care of such a region as illustrated in Fig. 1 in the 
following way. First we define the quantities [;; and A,; . 


‘ \1 if 72j is an internal point ah 
T,; = , (<9) 
0 if 2j isnot an internal point) 
J 
A; =1-—T;;. (26) 
Our computational formulas are written 
oe ae as (27) 
w 1 =a i _ wi; ‘Di + ; (28) 
wig = (Bi wi bisa (29) 
u Bio = 0; Uys — Be — 8. (30) 
B +1 1 
b . co | re 1 (31) 
w' 
RB 1; _ wi! 1, +1 
b i ve _— | Ag 9 (32) 
wr 
Bit? — wird 
at ‘oe! TTC (33) 
w; 
te Oe Meee Picasa (34) 
8 ee We Beciies (35) 
ek, =e... Bue #6 hx 
+ Bi" ‘6, stay Qigua ot Bj; s 1 Ajje1 (36) 
‘ ) 


i Osis Baie FB Rise Bec 
ae te re SO De ee ee Be 


where S;; is a term arising from sources and the values of 6 accompanying non-zero A 


are boundary values. 


E = D bia {(Bi;" a Bi; **") OF oil’, 1 1 (37) 
o4 


ij 


where 6* denotes a guessed value. 


E..-u I — wi,’ i — wis r 
si ad ee ' g df 1 Gi—13—~-1 5 ¢-14-1 (38) 
wi 
0, = Ji — ined Freee Weve — bi Osa. T iss; — bi" OsssisrT sais - (39) 


Some of the [’s and A’s are redundant here since the computed factors would be zero 
anyway. It is understood that the formulas are to be applied only when 7j is an interior 
point of the region. 

The calculational procedure for each iteration is as follows. We set our program 
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Fig. 1. A typical mesh. 


to run through the whole mesh as shown in Fig. 1 from the upper left hand corner to 
the lower right hand corner, computing and storing the b’s and g’s for the entire mesh 
using Eqs. (27) to (38). Then we run back through the mesh from lower right to upper 
left computing the 6’s using Eq. (39). 

4. The convergence condition for linear problems. Let us define the errors ¢ and ¢*, 


O.atc = Pais + €, (40) 

_ ee om Desas + r. (41) 

Values of @,....ea are specified for the entire mesh and put in as the 6*’s of (37). Then 
the @..,. are obtained from (39) in the course of the computation. Substituting (40) 


and (41) into (19) we obtain 
Bei. '€; li~a “TT B} .” 1j +6," 1j+1 


+ 2 + Rs + i HO es 


(42) 
fe ot ghee ee = —" ~~ ee ic 
+ (2 = Ge eile 
We can write this equation symbolically as 
Cit = Died . (43) 
In matrix and vector notation, 
Ce De*. (44) 


The solution for ¢ is written symbolically 


‘ C"' De* (45) 
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or 

e= Ke*, (46) 
where 

K = C" D. (47) 


Therefore we see that the convergence condition can be written in terms of the norm 
|| K || of K. It is 
Tah <1. (48) 


Thus if we can compute || K ||, we can tell whether or not a given problem converges. 


Let us take as a definition of the norm, the following expression 


| K || = Max Es (49) 


‘ i 


where the v; form a complete set of linearly independent vectors in the space of the 
problem. The process for obtaining Kv; from v; is identical to the process for obtaining 
e from e* by the use of (42). Equation (42) is similar to our original problem as given 
by (2); however, now we already have the 4’s on the left hand side of (42) so that it 
automatically satisfies the consistency conditions (16) and (17). Therefore, no iterations 
are required here. 

We thus investigate the convergence of a particular problem in the following way. 
First, we select the set of vectors v; . Then using (42), we compute the vectors Kv; . 
Then using Eq. (49) we compute the norm || K ||. If the convergence condition is sat- 
isfied, we proceed to obtain the solution. 

5. Numerical examples. As our linear example we will take the heat flow problem 
discussed by Baker and Oliphant [3]. The differential equation is 


; = & 


76 . 06 at 
oot aa, (50) 
Ox OY ot 


where a is a constant. The left hand side can be differenced either 


' So: mene ia1 — 40;; - 
1.6 = Sut Oisy + oF Baie — She (51) 
(Aa 
or 
L.0 = a Be ae us — 0 / (92) 


For our actual scheme, we take a weighted superposition of the schemes illustrated 
in (51) and (52). Thus, 
2 06 


<i (53) 
ol 


a,L,.6+a,L,6 =a 


where 
a, ta, = 1. (54) 


We now make a direct comparison with the linear problem worked by Baker and 
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Oliphant. For simplicity, we set (Ar) = (Ay). The time derivative is replaced by a 
three-level formula [4] 


00 363; — 4677" + 675" 


aS” 2(Ai) (55) 
Substituting (55) into (53), we obtain 
n—2 m—1) 
a,L,6+a,L,0 + B@=a” TT ha ; (56) 
where 
3 7 
P= ~ aah) (57) 


A simple algebraic manipulation shows that the Baker-Oliphant scheme for (Ar) = (Ay) 
can be obtained by setting 


es = Bian)? ? (58) 
a = Ban? (59) 


We observe that (54) is satisfied by (58) and (59). 
As in the particular example considered by Baker and Oliphant, we choose 


a’ (At)/(Ax)* = 1.02392228 so that the amplitude is diminished by a factor 10~'’’® 
at each time step. For this example, letting (Ar) = 1/12, we have 
B = —210.953511, (60) 
a, = 2.3652297, (61) 
a, = —1.3652297. (62) 


We carried this out on an IBM 704. The calculation gave the same result as the Baker- 
Oliphant calculation. It was observed by them that the solution for a single time step 
is exact and requires no iterations. Equivalently, we observe here that the 6};"'*' and 
8'*"'-* of Eqs. (34) and (35) are precisely equal to the B‘;*’*’ and Bi}''~’, respectively, 
so that the coefficients of 6*_,;,, and 6*,,;_, in (37) are precisely zero. 

Thus, we see how in this case our method reduces to the Baker-Oliphant method 
as a special case. We can look at the relationship between the two methods in the fol- 
lowing way. We obtain more generality by introducing the iterations on each time 
step. When we specialize to the Baker-Oliphant method we have no need for the iterations 
because the solution is obtained correctly in one step. 

To illustrate how our method is applied to a scheme which does not permit the Baker- 
Oliphant factorization, we considered an example with weights 


a, = 2/3, (63) 
i: = Ee: (64) 


As a criterion for iteration we required that successive iterates be within one part in 10° 
of each other. We found that for this example four iterations were required per time 
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step. In situations in which we can obtain a very good extrapolated guess for the suc- 
ceeding time step, we expect not to have an excessively large number of iterations, 
for example, some 5 or 6. This expectation has been borne out in a number of calculations 
of the above type. However, even in differencing schemes to which the Baker-Oliphant 
factorization does not apply precisely, it is advisable to choose weights closely cor- 
responding to their method in order to reduce the number of iterations as much as 
possible. In calculating linear examples, the program runs at about 13 milliseconds 
per cycle point. 

For our non-linear example we considered the same example 
and Oliphant, namely radiation flow in a material medium. The equation is 


as that used by Baker 


a, _ 060 
V'(6°) = 16K —. (65) 
ot 


With the present method, we need not make the change of variable required in the Baker- 
Oliphant scheme, but we can simply linearize the second partial derivative term. We 
use a linearization which is similar to that used by Bruce, Peaceman, Rachford, and 
Rice [2] in their one-dimensional examples. The resulting difference scheme can be 
written: 
(Bip Os-15-1 + Bip Os-ri4 HF Bi Ossian FH Bip" Oss 1541) 
+ a,(Biz 0:15 + Bij 0:41; + Big Oi5-1 + Bij "8i541) 66) 
(00 


. : (67-7 — 467°" 
— [(a,Bii(2) + a,B;i(+) + 8)]0;; = 16K a” im ; 


where 
Br” = f(0* , 0*,) = (0*)* + (0*)(0*,) + (0*)(0*,)” + (6%,)°, (67) 


except fori = m, j = n, 


Byi(x) = (O45 , OF ass) + AO% , OR asas) + HOF , OB asa) + HOF , OF s+), (68) 


i 


Bii(+) = (05 , 0F1;) + (O48 , OF15) + (05 , OF -1) + f(O4% , 6% 41), (69) 


and 
i) a 
aah (70) 


ee Tan 


In the non-linear case we see that the 6*’s occur on the left-hand side as well as the 
right-hand side. The procedure to be followed in accordance with the analysis of Sec. 4 
would be to put a guess 6** into the coefficients given by Eqs. (67) to (69) and iterate 
6* in Eq. (37) down as in a linear equation, and then to make another guess 0** and 
so on. In practice we find it much better to guess 6* for both (67) to (69) and (37) at 
the same time. The latter procedure was followed in the example discussed below. 

For our sample calculation we picked 8 Ax” = B Ay’ = —1.5 and used 11 X 11 
interior mesh points. For initial conditions we chose 


—_ 9 
Ax, y) = 1 —sin (=) sin (2=2), (71) 








1961] TWO-DIMENSIONAL TIME-DEPENDENT DIFFUSION PROBLEMS 229 


and we maintained the boundaries at unit temperature. This is the identical example 
chosen by Baker and Oliphant. Initially some 6 or 8 iterations were required per time 
step, but after 16 time steps, only about 2 iterations per time step were required for 
three-place accuracy. We made the initial guess at each time step by a linear extrapola- 
tion from the previous two times. The program runs at about 47 milliseconds per cycle- 
point. 

As with the Baker-Oliphant method, care should be exercised to prevent the guessing 
of a negative temperature, as the occurrence of temperatures of different signs produces 
an instability which causes the solution to diverge. 

As was observed in the above calculations, the convergence of the iterations was 
rather good. We can discuss why this should be so in the linear case on the basis of the 
analysis of Sec. 4. We computed the norm of (49) for a number of sample matrices. 
It was found that the convergence is best for the heaviest weighting of the diagonal 
element of BY; . On the left-hand side of Eq. (56), the first two terms weight the diagonal 
and off-diagonal terms equally. In steady state problems, these are the only terms 
present and for such problems we get fairly good convergence. For the time-dependent 
problems, the third term adds on in such a way as to increase the weighting of the 
diagonal elements and thus to enhance the convergence greatly. Also, in time-dependent 
problems, taking extrapolated guesses helps considerably. 
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Lectures on differential and integral equations. By Kosaku Yosida. Interscience Publishers, 
Inc., New York and London, 1960. ix + 220 pp. $7.00 


This text is the English edition of a volume published originally in Japanese. It is intended to be a 
self-contained exposition of the theory of ordinary differential equations and of integral equations. 
The subject of partial differential equations lies outside the scope of the book. Chapters, 1, 2 and 5 
deal with ordinary differential equations, while Chapters 3, 4 and 6 are concerned with integral equa- 
tions. Chapter 1 begins by treating the initial value problem dy/dx = f(z, y), y(zo) = yo, for a single 
equation and for a finite system, both in the real and complex domain. After dealing with linear differ- 
ential equations of the nth order, the chapter ends with a succint treatment of second order differential 
equations of the Fuchs type, with the special cases of Gauss, Legendre, and Bessel being analyzed in 
detail. Chapter 2 covers the Sturm-Liouville boundary value problems for linear second order differential 
equations, Green’s functions; generalized Green’s functions, the Hilbert-Schmidt theory of integral 
equations with a symmetric kernel, and Liouville’s method for obtaining asymptotic expressions 
for the eigenvalues and the eigenfunctions. Chapter 5 is devoted to an elementary exposition of the 
general expansion theorem (Weyl-Stone-Titchmarsh-Kodaira) relative to the differential equation 
y’ + (A — q(xz)) y = 0,a < x <b, where g(z) is a real valued continuous function and no assumption 
is made about the behavior of g(z) as z — a or as x - b (this is the general “singular” case). Let y;(z, ), 
yo(x, ) be a fundamental system of solutions determined by initial conditions at a number c with 
a<c <b: y(c, A) = 1, yilc, ) = 0; yo(c, A) = 0, yc, X) = 1. Then, for “appropriate boundary 


conditions at x = a and at x = ,b’’, every real valued continuous function f(z) ina < x < b with 
Safx)dz < + © can be “expanded” as follows: 
+o 2 u b 
fa)= [| ah dD f ule,w draw | suls, w asp. 
—2 i,k=1 a 


This is the Weyl-Stone expansion theorem, which has been completed by Titchmarsh and Kodaira 
by giving an explicit formula for the ‘‘density function’”’ p;.(u). This general expansion theorem enables 
the author to give a unified treatment of classical expansions in terms of special functions, such as the 
Fourier series expansion, the Fourier integral, the Hermite polynomials expansion, the Laguerre poly- 
nomials expansion, and the Bessel functions expansion, which are worked out in detail as special cases 
of the general theorem. For simplicity, the casea = — ~,b = + © is considered in the proof of the 
general expansion theorem, which is obtained as a limiting case of the previously obtained Hilbert- 
Schmidt expansion theorem for the “‘regular’’ case of a finite interval. Chapter 3 is dedicated to the 
theory of Fredholm integral equations of the second kind and Chapter 4 to the theory of Volterra 
integral equations of both the first and the second kind. Chapter 6 deals briefly with non-linear Fredholm 
and Volterra equations of the second kind. The author has succeeded in achieving an amazingly com- 
pact presentation of the material. This book is remarkable both for its conciseness and for its readability. 


J. B. Draz 
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ON ‘“‘TRANSCRITICAL” AND ‘‘SHYPERCRITICAL” FLOWS 
IN MAGNETOGASDYNAMICS* 


BY 
R. SEEBASS 
(Cornell University) 


Introduction. A category of magnetogasdynamic flows that has attracted con- 
siderable interest involves the steady motion of a perfectly conducting, inviscid, com- 
pressible fluid with magnetic field everywhere aligned with the flow direction. This 
type of flow must result if the magnetic field and flow direction are aligned anywhere in 
plane and axially symmetric configurations. It is clear that such flows are isentropic, or 
in the case of flow from a uniform state, homentropic. 

The system of equations governing motion in this category may be either elliptic 
or hyperbolic, depending upon the local values of the Alfvén number, A, and the Mach 
number, M. The Alfvén number is defined here as the ratio of the flow speed to the 
propagation speed of Alfvén waves, while the Mach number is the usual ratio of flow 
speed to the speed of sound in the absence of magnetogasdynamic effects. These elliptic 
and hyperbolic regions were first exhibited by Taniuti [1] in a diagram similar to Fig. 1, 
and later by Kogan [2] and Resler (see Sears [3]). There are two regions where the equa- 
tions are hyperbolic, one of which is subsonic, sub-Alfvénic (A < 1), the other being 
supersonic, super-Alfvénic (A > 1). In these regions the proper family of characteristics 
must be chosen. Sears [4] has pointed out that for the subsonic case waves which propa- 
gate upstream are the correct choice, while in the supersonic region it is the usual down- 
stream-propagat ing waves. 

Along PR, where A® + M’ = 1, the propagation speed of small magnetosonic dis- 
turbances vanishes. The elliptic region OPR where the velocity is less than the propaga- 
tion speed will be referred to as subcritical. The subsonic, hyperbolic region PRQ where 
the velocity is greater than this speed will be called supercritical. Near PR the flow will 
be termed hypercritical. Note that both transonic and trans-Alfvénic flows are, under 
this nomenclature, transcritical. 

McCune and Resler [5] have derived the linear theories for this type of flow. Just 
as the results of ordinary linearized gasdynamics become invalid in the transonic regimes, 
McCune and Resler’s results break down in the transonic, trans-Alfvénic and hyper- 
critical regimes. The object here is to study the flow in these regimes where the motion 
is fundamentally non-linear and the non-linear terms of the equations must be retained. 
Of particular interest is the hypercritical regime where elliptic flow joins with hyperbolic 
flow having forward-facing characteristics. 

It is found that for the flow under consideration the hodograph transformation 
can easily be effected, the result being two second-order, linear, partial differential 
equations analogous to Chaplygin’s equation and the potential equation for ordinary 
compressible flow. Two elementary solutions of these equations, corresponding to source 
and vortex motion, are discussed. One displays smooth transition through the hyper- 
critical and trans-Alfvénic regimes. Separation of variables is possible, and the resulting 
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ordinary differential equations are discussed briefly. In the special case where the ratio 
of specific heats, y, is 2, the equations are of the Fuchsian type with five singularities. 

Equations of motion. The equations governing the steady flow of an ideal conductor 
in the presence of a steady magnetic field are 


div pq = 0, (1) 

grad qg°/2 — qX curl q + (u/4rp)HX curl H + (grad p)/p = 0, (2) 
div H = 0, (3) 

curl (q X H) = 0, (4) 

q-grad S = 0, (5) 


where p, p, S, » are the pressure, density, entropy and magnetic permeability, and q 
and H the velocity and magnetic field vectors. In plane flow, Eq. (4) requires that 
q X H be constant; thus if q and H are parallel anywhere (e.g., at infinity), they are 
parallel everywhere. This can be expressed as 

pq = fH, (6) 
where f is some scalar function. Equations (1) and (3) then give f = f(W), where y is 
the stream function. Substituting Eq. (6) into (2), forming the scalar product with 
q, and noting Eq. (5) gives Bernoulli's equation 


q/2+ | dp/p = constant ( 


along a streamline. Assuming uniform conditions at infinity, f becomes constant, the 
flow is homentropic, and Eq. (7) holds through the flow. Combining the gradient of 
Eq. (7) with Eq. (2) and (6) then yields 
curl (1 — A~*)q = 0, (8) 
where A denotes the Alfvén number, (42p/u)*¢/H. 
The characteristics for this system, first given by Taniuti [1], can be expressed in 
the form 
a[A° + M* — 1)[(dx)’ + (dy)*] = A®*[u dy — v dx]? (9) 
in the physical plane, and 
(M* — 1)[A’? + M?’ — 1][udv + v du)’ = (A? — 1)[u dv — v dup? (10) 
in the hodograph plane, where a* denotes dp/dp and M, q/a. From Eq. (10) the diagram 
in Fig. 1 can easily be deduced. Equation (9) shows that as A~ + M* — 1, the character- 
istics become parallel to the streamlines, and as A — 1 or M — 1, they become per- 
pendicular to the streamlines. 
Hodograph transformation. Equation (8) suggests the introduction of a potential 
function, ¢, such that 
grad g = (1 — A’°)q. 


With this equation, the usual definition of the stream function, Bernoulli’s Equation (7), 
and the continuity Equation (1), the transformation to the hodograph plane can be 
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effected in a manner identical to that for the analogous equations in ordinary compressible 
flow (see e.g. von Mises [6]). The first-order, linear, differential equations 





i a A (11) 
66 (A> + M* — 1) p dq’ 

dp _ ce — 1)(A* — 1) dy (12) 
0q A* pq 06’ - 


are the result, subject to the condition that the Jacobian of the mapping does not vanish. 
In terms of y alone this requires that 


a(x, y) I dy\’ (A* — 1)q_ (dp\’ ; 
‘ = M’? »( ) - ng a 0. 13 
0(q, 6) pq. [ar - 06 A? ry Mw oq (13) 
Limit lines, i.e. curves on which J = 0, can occur only in the two hyperbolic regions. 
In contradistinction to ordinary compressible flow, the Jacobian 
dy, W) » (A? — 1) 


= a ° I 
O0(q, 8) pq A’ 








can be zero even though J is not zero. This is a result of the definition of ¢ for this flow. 
Elimination of the potential function from Eqs. (11) and (12) gives, for a perfect gas, 


rat 4. * ns SPS ’ sail ’ 7? 
iu? — 1) 4+" oad ~(2 -1(42+M?- 5a 
. (14) 
ail 1 dy 
— {(1 + M’)(A? — 1)? — MTA )+1-3a} 2 0, 


which reduces to Chaplygin’s equation in the limit of no magnetic field, A —~ ©. The 
equation for ¢ alone is 





.2[A?> + M* — 1) 0a 9” 
a — m4 f Nee. gy — ya? — 1) <4 
q 00 oq (15) 
4\ 2 2 2 l dg 
+ [1 + yM)\(AL — 1) + MM — 1) - de = 0. 
To return to the physical plane the following coordinate relations are required 
ens | fo  emtly ~ 2 8 ay, (16) 
qLA°—-1 
1 | cos A’ ’ ” 
dy = cm + — a7 -15" Ode}. (17) 
d 


Elementary solutions. [Equations (14) and (15) have the “elementary”’ solutions 
vY = Kéand¢g = K@, where K isa constant. The corresponding solutions in the physical 
plane could have been obtained directly from Eqs. (1) and (8), without recourse to the 
hodograph transformation. 

The first of these solutions, which might be termed “magnetic-source” flow, has 
a two-fold mapping to the physical plane and the limit line 17 = 1. Along this line the 
two solutions, one subsonic and the other supersonic, join. These are depicted in Fig. 1 
by the curves Ia and Ib. The arrows indicate the direction of increasing radius. Equa- 
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Fig. 1. Diagram of elliptic and hyperbolic regions of flow. 


tions (16) and (17) give the coordinate functions x and y, which are identical to the 
non-magnetic results and will not be given here (see e.g. [4]). The subsonic flow, Ia, which 
undergoes smooth transition through the hypercritical and trans-Alfvénic regimes, 
is shown in Fig. 2. 
The second solution, a ‘‘magnetic-vortex’”’ flow, has the limit line A? + M’ = 1 and 

the coordinate function 

in geal a2 _ , KA" 

r=("e+y)" = +A? oh 
In this case the mapping to the physical plane is threefold, the flows being indicated 
in Fig. 1 by the Curves IIa, IIb, and IIc. One set of flows joins at the limit line, the 
other at a branch line 1/J = 0. The flow corresponding to IIb is shown in Fig. 3. 


The special case y = 2. Assuming a solution to Eq. (14) of the form y = [(@)F(r) 
results in 


f =ce" + ce (18) 


and a second-order, linear, differential equation for F, where r = q°/qi,; and k is the 
separation constant. The coefficients of this equation are, in general, complicated alge- 
braic functions of 7. With y = (n + 1)/n, n a positive integer, the coefficients become 
polynomials in 7, the degree of which depends on n. In the simplest case n is 1 and the 
equation for F becomes 

ral r(r — 8)’ ar 3k? (r — B)(1 — 37) p _ 


dr 4 ar 1)’ 0, (19) 


(r — B/3)(r — 1) dr 
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Fic. 2. Streamlines and characteristics for ‘‘magnetic-source’’ flow. 


where 1 — 6 A* evaluated at r = 0. Here there are five singularities, all regular, 
and thus the equation is of the Fuchsian type. In terms of the extended Riemann 
P-function [7], 


|. (20) 
\+hk/2 2 —1 2 +3h/2 


The singularities 8/3 and 8 correspond to A* + M* = 1 and A® = 1 respectively. Results 
analogous to Eq. (18), (19) and (20) may also be obtained for the potential function. 
They are more complicated in that the exponents at two of the singularities are func- 
tions of @. 

An interesting conclusion can be drawn in this case with regard to the flow behavior 
as A~X + M*° — 1. The singularity at 7 = 8/3 can be shown to be apparent, and therefore 
both solutions are analytic about this point. This insures further that J will be bounded 
and generally non-zero. In the hypercritical regime then, for y = 2, the flow behavior 
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Fic. 3. Streamlines and characteristics for “‘magnetic-vortex’’ flow. 


must be smooth. On physical grounds it might be conjectured that this will hold for 
other values of y. In contrast, 0 and 1 are the only real values of k for which r = 1 is 
an apparent singularity. 


It is possible to solve Eq. (19) in terms of rational functions of 7 for particular values 


of k. The procedure is too lengthy for presentation here. For k = 1, the two solutions are 
; l 
F = —7y 
ro wr = 2 
and 
r*(27? — (8 + 3)r + 28) 


7-6 


Solutions of this class should provide insight into the flow behavior in the hypercritical 
regimes analogous to that obtained from the solutions of Ringleb [8] or Temple and 
Yarwood (see [9]) in ordinary gasdynamics. 
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In the limit 8 — 0 Eq. (20) reduces to the Riemann P-equation with the solution 
PF  iiinaaiala , 


=i => Bk’) —1 + 3k 3k°)'? et 
R( | — 3k + (1 + 3k! 1 + 3b + (1 + 3¥/) + (1+ 38"; *) 


_ — 


where F is the hypergeometric function. For k = 1 a solution is 
3/2: ~1/2 -1 
y=r7 sin @, g=7r (1 — 7)" cos 0. 
Since 8 = 0 the motion is entirely super-Alfvénic. The limit line and streamlines are 


similar to those of Ringleb flow. 

Conclusion. Through this investigation the existence of smooth flows in the 
transonic, trans-Alfvénic, and hypercritical regimes has been determined. The flow 
behavior in the hypercritical regime will generally be regular for the special case y = 2. 
The new and interesting phenomena associated with the flow of a highly conducting 
gas in the presence of a strong aligned magnetic field have been exemplified in two 
simple solutions. 

Acknowledgment. The author is indebted to Professor W. R. Sears for suggesting 
this investigation and for his guidance through its course. 
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BOOK REVIEWS 


(Continued from p. 230) 


Mathematical methods for digital computers. Edited by Anthony Ralston and Herbert 
S. Wilf. John Wiley & Sons, Inc., New York, London, 1960. xi + 293 pp. $9.00. 


This collection of articles constitutes a handbook for digital computer programmers, giving al- 
gorithms for the solution of standard problems which a programmer should be able to translate into 
any desired code. Each article adheres to a consistent format, giving the purpose of the program, the 
mathematical background, the calculation procedure, the flow-chart and its description, specification 
of required subroutines, a sample problem, and estimates of memory and running time requirements. 
There is, naturally, much variation in the treatment of each of the above topics by the contributors, 
some giving detailed, some only large scale flow-charts, some more and some less mathematical back- 
ground. One might wish that the editors had requested from each contributor a FORTRAN or ALGOL 
program to accompany the flow charts and thus imposed a uniform standard in this respect. But even 
without this, the book will be greatly welcomed as a pioneering example of a type of work in great 
demand amongst applied mathematicians with computer orientation. 

There are 26 articles, divided into 6 groups. The first article, alone in its group, deals with the 
generation of elementary functions. In it, E. G. Kogbetliantz summarises his contributions previously 
published in the IBM Journal of Research and Development; he discusses power series and rational 
approximation in a somewhat concise and difficult to read paper which, however, merits the attention 
it demands. 

The second group, on matrices and linear equations, consists of six papers. Direct methods in matrix 
inversion and related topics are discussed by Alex Orden. He gives an omnibus program which (1) finds 
the inverse if it exists, (2) solves systems of linear equations, (3) finds the value of a determinant, 
(4) determines rank and (5) locates a nonsingular minor. The matrix has, apparently, to be small enough 
for each of its elements to be stored in memory, and dense enough for this to be reasonable—there is no 
discussion of tape-techniques, a criticism applicable to most of the matrix papers in the volume. In the 
second paper, R. Van Norton describes a Gauss-Seidel program without special acceleration devices; 
background theory and convergence criteria are given. Next, F. S. Beckman presents the conjugate 
gradient method; this iterative process which, unlike Gauss-Seidel, gives the solution after a finite 
number of steps, was invented by Hestenes and Stiefel in 1953 and in his careful exposition of the 
algorithm the author shows its accuracy and speed to be high when applied to sparse matrices, but not 
otherwise. Herbert S. Wilf first published his method of rank annihilation for matrix inversion in the 
SIAM Journal in 1959; the next paper, by him, is devoted to this method which depends on a formula 
giving the change of inverse in terms of the change in value of one element. The von Neumann and 
Ulam Monte Carlo method for matrix inversion is expounded by Florence Jeanne Oswald, with a sampie 
problem. The last paper in this group is by John Greenstadt, on the Jacobi method for the determination 
of the eigenvalues of a matrix. The recommended version is the ‘threshold’? method of Pope and 
Tompkins (J. Assoc. Comp. Mach., 1957). The presentation is confined to symmetric matrices, which 
is a pity since the author and others have recently made progress in the difficult general case; perhaps 
it is not yet possible to present a general computer algorithm for general non-symmetric matrices— 
it should be possible, however, if they have a complete set of eigenvectors and can thus be diagonalized. 

The next group of four papers is on the numerical solution of ordinary differe.tial equations. First 
of all, Anthony Ralston gives an informal discussion of multistep methods for initial value problems of 
first order, including truncation error analysis, variable step-length procedures and stability considera- 
tions. Next, Michael J. Romanelli describes Runge-Kutta methods for first order equations, both with 
and without the Gill modification. To exemplify the numerical solution of boundary-value problems, 
Eugene L. Wachspress chooses a Sturm-Liouville problem; the tridiagonal system of linear equations 
arising from the finite difference scheme is solved by Gaussian elimination; round-off error propagation 
is discussed, convergence and trunction errors are not. Next, J. Certaine is concerned with the solution 
of ordinary differential equations with large time constants; these are well known to present numerical 


(Continued on p. 252) 
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A DUALITY THEOREM FOR NON-LINEAR PROGRAMMING* 


BY 
PHILIP WOLFE 
The RAND Corporation 


Summary. A dual problem is formulated for the mathematical programming 
problem of minimizing a convex function under convex constraints which reduces to 
the classical dual problem in the case of linear programming problems. Duality theorems 
are proved regarding the relationship between the problem and its dual. 

1. Introduction. A duality theorem in mathematical programming is, generally 
speaking, the statement of a relationship of a certain kind between two mathematical 
programming problems. The relationship commonly has three aspects: (a) one problem— 
the ‘‘primal’’—is a constrained minimization problem, and the other—the “dual’”—is 
a constrained maximization problem; (b) the existence of a solution to one of these 
problems ensures the existence of a solution to the other, in which case their respective 
extreme values are equal; and (c) if the constraints of one problem are consistent, while 
those of the other are not, there is a sequence of points satisfying the constraints of the 
first on which its objective function tends to infinity. 

Dennis [1] and Dorn [2, 3] have formulated dual problems for the primal problem 
of minimizing a convex function f of several real variables under linear constraints, 
and shown (b) above in the case that f is quadratic, or just differentiable but strictly 
convex. This note formulates a dual to the problem of minimizing a convex differentiable 
function under non-linear, concave constraints; this dual problem reduces to theirs in 
the case of linear constraints. 

It has not been possible to establish (b) and (ce) in full, and there is evidence that 
they do not hold as stated. Theorem 2 below derives the existence of a solution to the 
dual from that for the primal, and Theorem 3 establishes (c), in part, in the case of 
linear constraints. Thus the problem studied here does not seem to enjoy the complete 
symmetry of the notion of duality which is found in the linear case. 

The known state of affairs regarding duality for this sort of problem is summarized 
in the table of Sec. 4. Several questions remain open. The utility of all this is discussed 
in Sec. 5. 

2. The problems. Let f be a convex, differentiable function of x = (x, , «++ , 2,); 
for each 7 = 1, --- , m, let g;(x) be a concave, differentiable function. For any function 
¢, let V@ denote the gradient of ¢, 


(22. at a) 
oz,’ ? 0z,J- 


Primal problem Minimize f(x) subject to 





g(x) > 0, i= 1, eS 5 aoe (1) 
Dual problem Maximize f(x) — > u.g,(x) subject to (2) 
Vif(x) = p 2 u;V g(x), u> 0. (3) 


*Received August 25, 1960; revised manuscript received February 20, 1961. 
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A solution of one of these problems will be a point, x or (x, u) respectively, which 
achieves the extreme value sought. The problems will not in general have solutions; 
about the only relation we can state in their absence is that of Theorem 1 below. 


Theorem 1. Let V be the infimum of f(x) under the constraints (1), and v be the 
supremum of (2) under the constraints (3). Then 


v<V. 
Proof. Let x* satisfy g;(x*) > 0 (all 7), and (2, u) satisfy the constraints (3). We 
then have the chain of relations 
f(x*) — f(z) > 3=Vf(x)(2* — 2) 
= p u;Vg(x)(x* — 2) 


IV 


> ulgie*) — g.(x)] 


— 7 U:gi(X); 


the first is due to the convexity of f; the equality, to the constraints (3); the next in- 
equality to the concavity of the g; ; and the last inequality to the non-negativity of u, 
and g;(x*). It follows that 


IV 


f(x*) > f(x) - 2. U.gi(x), 


which proves the theorem in the case that both sets of constraints are consistent. Allow- 
ing the convention V = + © in case (1) are inconsistent, and v = — © in case (3) 
are, the theorem follows. 

3. The duality theorems. The “constraint qualification” of Kuhn and Tucker 


{6, p. 483] will be assumed from now on. Stated for the primal problem at hand, it runs 


Let x satisfy (1), and let dx be any vector differential such that Vg;(x) dx > 0 for 
all 7 such that g;(z) = 0; then dz is tangent to some arc contained in the set of all 
x satisfying (1). 
The qualification serves to rule out certain singularities which might otherwise occur 
on the boundary of the constraint set. The constraints g; will meet the qualification if, 
for example, the constraint set has an interior point relative to the convex determined 
by those constraints which are linear; and thus a set of entirely linear constraints meets it. 
Under the constraint qualification and the properties that have been assumed for f 
and g; , the Equivalence Theorem of Kuhn and Tucker [6, Theorem 3] reads 


rr . 0O- . . . “oe “fe 0 
rhe point x’ is a solution of the primal problem if, and only if, x” and some w’ con- 
stitute a saddle point for the Lagrangian function 


L(x, u) = f(x) - » Se U;:gi{X), (4) 


that is, 


for all w > O and all x. 
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Our principal duality theorem is 


Theorem 2. If x° solves the primal problem, then there exists u° so that (2°, u°) 
solves the dual problem, and the extrema are equal.* 

Proof. The function L(x, u) is convex in x for any u > 0. It consequently satisfies 
the following relation, which asserts that its graph as a function of x lies above any of 


its tangent planes: 
L(y, u) — L(x, u) > V.L(2, wy — 2). 

Now V.L(z, u) = 0 for any (x, wu) satisfying (3), so that if (x', u) and (x’, u) both 
satisfy (3) we have both L(2', u) — L(x’, u) > 0 and the reverse inequality, whence 
L(x’, u) = L(x’, u); in other words, L(x, u) is independent of x for (x, u) satisfying (3). 
Consequently. 

L(x’, u’) = Max {L(2°, u) | u > 0} 


> Max {L(x’, u) | (x", u) satisfy (3)} 
= Max {L(z,u) | (2, u) satisfy (3)} 
> Liz’, w), 


so that L(2°, u°) is the maximal objective value for the dual problem. This proves the 
first part of the theorem. 
Finally 


L(a’, u°’) = f(a’) — >> ulg,(2’) = f(z’), 


because each u,g,(2°) = 0; if this last statement did not hold we would have u; > 0 
and g,(x°) > 0 for some i, and then L(x’, u°) could be increased by decreasing u? , in 
contradiction to the saddle-point property (5). The proof is complete. 


Theorem 3. If the primal problem has only linear constraints and these are in- 
consistent, and the constraints of the dual problem are consistent, then the supremum 
of the objective function (2) of the dual problem is + @. 

Proof. Let (x, u) satisfy the constraints (3) of the dual problem, and let 
g,(2) = A‘x — 6, (all 2). Since these latter constraints are inconsistent, there exists 
(4, Lemmas 4 and 5] u* > 0 such that u*A‘ = 0 (all 2) and > u,b; > 0. It follows that 
the point (z, wu + tu*) satisfies the constraints (3) for all ¢, and that f(z) — 
} (u; + tu*)g,(x) ~ + ~ ast + o. 

4. The state of affairs. The following four possibilities exist for both the primal 
and the dual problem: (1) the problem has a solution; (2) the constraints of the problem 
are consistent and the function to be extremized is bounded in the direction of extremi- 
zation, but the problem has no solution—the bound is not attained; (3) the constraints 
are consistent but the function is not bounded in the direction of extremization; (4) 
the constraints are inconsistent. 

The table below summarizes the known results in terms of the possibility of a given 
pair of outcomes for the two problems. Theorems 1 and 2 account for the impossible 
outcomes; examples exist for the outcomes known to be possible. Three cases are un- 


*I am indebted to W. 8. Dorn and R. E. Gomory for pointing out an error in an earlier version 


of this theorem. 
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settled; an example of a problem whose dual was bounded but had no solution would 
help clear them up. With the exception of the second row and the second column, which 
would constitute impossible categories for linear programs, and the annoying example 
(6), the table is just the same as it would be if both problems were linear. 


TABLE. 


Outcomes of primal and dual problems. 


PRIMAL 


Has Bounded, Unbounded Inconsistent 


DUAL solution no solution 
Has YES ? NO YES 
solution (5 (1 (6, 3 
Bounded, NO ? NO ? 
no soln. 2 (1 (3 
Unbounded NO NO NO YES 
ae. (] (1 (5 
Inconsistent NO YES YES YES 
‘2 } 5 (5 


(1) Established by Theorem 1. 

(2) Established by Theorem 2. 

(3) “NO” if the constraints are linear, by Theorem 3. 

(4) Illustrated by the example: m = 1, n = 1, f(x) = e?, g(x) = 0. 

(5) Illustrated by readily constructed linear programming problems [4, pp. 57, 58). 

(6) Illustrated by the example: m = 1, n = 1, f(x) = 0, g(x) = —e?. 

5. Critique. Suppose that the constraints of the primal problem are linear 

g(x) = Ayx — Db; (¢ = 1, «++, m), 

where A, is a row of the m by n matrix A, and b, is the ith component of the vector b. 
The constraints may then be written 


Ax > b. 
The dual problem assumes the form 
Maximize f(z) — uAx + ub subject to 
Vi(«) = uA, u> 0. 
(Dorn [2] gives the objective function in the equivalent form 
f(z) — Vf(x)e + ub.) (6) 
In the case of quadratic programming with linear constraints, 


f(x) = px + $2Qz, 
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where Q is an n by nm symmetric matrix. The objective function for the dual problem 
becomes 


px + 3aQx — uAx + ub = —}2Qz + ub, 
and its constraints 
p+2Q = uA. 


Finally, of course, in the case of linear programming—quadratic programming with 
( = 0—we obtain the classical dual problem 


Maximize ub subject to 
p=uA, u > 0. 


A basic difference between the non-linear duality results and those for linear pro- 
gramming is that in the non-linear case the function f of the primal problem appears 
not only in the constraints of the dual, as expected, but remains involved in its objective 
function as well. This fact seems to diminish somewhat the usefulness of the dual prob- 
lem for computational purposes: the primal should commonly be easier to solve. (How- 
ever, the “‘cutting plane’ method [5] for non-linear problems, formulated in dual terms 
for computational purposes, turns out to be effectively an algorithm for the direct 
solution of the dual non-linear problem.) 

In certain cases the explicit appearance of x in the objective of the dual problem 


can be transformed away. If the mapping y = V f(x) is inversible (which will be the 
case if f is strictly convex), so that x h(y) may be written, then the dual problem 


in the ease of linear constraints assumes the form 
Maximize f(h(y)) — yh(y) + ub subject to 
y= uA, u> 0. 


Dennis [1] states the dual problem in this way, using the “Legendre transform”’ 
f(h(y)) — yh(y) of the function f. His formulation is more attractive when f if quadratic. 


The inversibility of f entails the non-singularity of Q. Letting R be the transpose of 


the inverse of Q, under the substitution of y for x the dual of the quadratic problem 


€ 


can be stated 
Maximize — }yRy + pRy + ub — 3pRp subject to 
y = uA, u> 0. 


It should be noticed that the “nice’’ properties regarding convexity that have been 
postulated for the primal are not inherited by the dual in the general case, or even in 
the case of linear constraints. The original form (2) for the objective of the dual problem 
is a convex function of x (one does not like to maximize such a function), while its form 
6) is neither convex nor concave; the set of (x, u) satisfying the dual constraints is not 
convex; and while the objective of the dual depends only upon u for points of this set, 
it is neither a convex nor a concave function of wu. 

Things are much better, however, in the quadratic-linear case. The constraints are 
linear and the objective is concave. The linearity of the constraints permits linear com- 
putations in handling the dual, and makes possible an entirely linear algorithm for 
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the solution of this problem. In the ‘‘simplex method for quadratic programming” [7], 
a point (x, u) satisfying the constraints of the dual is found for which x satisfies the 
constraints of the primal and such that the difference, u(Ax — b), between the objective 
functions of the two problems vanishes. By Theorem 1, such a point solves the problem. 
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A VARIATIONAL PROBLEM RELATING TO COOLING FINS 
WITH HEAT GENERATION* 


BY 
CHEN-YA LIU 
Carnegie Institute of Technology, Pittsburgh, Pa. 


1. Introduction. Cooling fins are used in heat exchange apparatus to increase the 
rate of heat transfer. To economize in the design, we wish to know the shape of the 
fin which gives the maximum dissipation of heat for a given weight of the fin. For pure 
conduction fins, a criterion for this optimum fin problem was proposed by Schmidt [1] 
and recently proved by Duffin [2]. Cooling fins are also used in atomic reactors where 
heat is produced inside the fin as a result of atomic reaction. The question immediately 
arises: What is the optimum fin gecmetry in such cases? The answer to this question 
becomes more important in view of the industrial trend of trying to develop airborne 
reactors where the weight limitation is the most significant problem. 

In the first part of this paper, the problem of cooling fins with heat generation is 
recast in a form suitable for treatment by the calculus of variations. The heat generation 
function which is not clearly known in our present state of knowledge is assumed as 
a function of the coordinate along the fin. The relation of the temperature to the heat 
generation is assumed to be linear. Euler equations are obtained by formal variational 
methods. Contrary to the case for pure conduction fins, the equations are not linear. 
General solutions to the Euler equations cannot be obtained in explicit form. However, 
sufficient conditions are derived for solving this optimum fin problem. 

The second part of this paper concerns the solution when the heat generation function 
is linearly dependent on the temperature only. The temperature is found to be a hy- 
perbolic sine function. This result is used to derive explicit expressions for the fin shape 
and for the maximum heat dissipation, for the cases of a rectangular fin and a circular fin. 

2. The variational problem. Under the assumption that the heat generated is 
linearly proportional to the temperature u, the governing equation for a cooling fin 


can be written as 


d --. Gu ’ 
- , = [ (7) — " r) | 
ym Ee au [p(x) — q(x)y(x)]u, (1) 


where the known functions p(x) and q(x) representing the surface convection and heat 
generation effects respectively are positive and continuous. 
The function y(x) related to the fin shape satisfies the conditions 


y(x) is differentiable (2) 
yx) >0O for O<a<b (3) 

and 
| y(x) de = K, a given constant (fin volume), (4) 


70 


where b is the fin length. 


*Received January 26, 1961. 
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Also u is to satisfy the two boundary conditions 


or 
~~ 


u=1 at z=0 ( 


) du 


y(x) a? 0 at t= bd. (6) 


Under these conditions, we seek two functions u(x) and y(x) and the constant number b 
such that the cooling effect of the fin H defined as 


H = [ [p(x) — g(x)y(x)]u(x) dx (7) 
70 
gives a maximum. 
Further we restrict the function y(x) such that 


p(x) — g(x)y(x) > 0. (8) 


Thus the cooling effect of the fin is ensured. 
The solution of the problem can be made to depend on the following integral 


E = i | (a4) +o = aa’ | dx. (9) 


Integration by parts and use of (1), (5) and (6) gives 


} 


> du? re i du) - 


70 /0 J0 


Substituting this in (9) shows that 
ab 
E= | (p — qy)u dz. 10) 
0 


In other words, if u satisfies (1), (5) and (6), then EF = H. 

Equation (10) suggests that we minimize the integral 7, treating u and y as inde- 
pendent variables, not supposing that wu and y necessarily satisfy Eqs. (1), (5) and (6). 
Note that under the conditions (3) and (8), the integrand in (9) is never negative and 
hence E has a finite lower bound. 

We proceed formally, considering a variation of E(u, y) due to a variation of u and y 
from the assumed minimizing values. First only u is varied. The Euler equation resulting 
for any admissible variation of u is the governing equation (1). Now consider a variation 
of E resulting from a variation in y. Thus 


b 29 
, du\* 2 
bE = [ ay| (2) = ai | dx = 0. 


Because of condition (4), we also have 
h 


6K = [ by dx = 0. 


v 
It follows that the integrands of 6 and 6K must be proportional. Thus 


(du/dx)? — qv’ = x. (11) 





1961] COOLING FINS WITH HEAT GENERATION 247 


Integrating Eq. (1) and use of (6) gives 


du _ F 9 
i “i -/ (p — qy)u dz. (12) 


In view of relations (3), (8) and (12), Eq. (11) can be written as 


du/dx = —(qu’ +r)”, (13) 


where \ is a constant (Lagrange multiplier). The minimum of £ is obtained from the 
functions u and y which simultaneously satisfy Eqs. (1), (5), (6) and (13). 

Because the non-linearity of Eq. (13) prevents us from finding explicit solutions, a 
general procedure for solving the optimum fin problem will be given before specific 
examples are discussed. 

CasE (a). Optimum fin. After solving Eqs. (1), (5), (6) and (13), we will obtain 
a one parameter family of solutions u(x, A) and y(z, b, A) which, if substituted in Eqs. 
(4) and (7), give K = K(b, \) and H = H(b, d). Hence we have H as a function of b 
and \ which in turn are connected by K. Maximizing H while holding K constant will 
give the optimum fin length. The optimizing parameter \ can then be evaluated from 
K(b, Xd). 

An alternate method which gives rather simple conditions for the optimum solutions 
is obtained as follow. In deriving the above variational result, the fin length b was 
held fixed. However, in the problem stated, b should be allowed to vary also. For this 
purpose we consider the function F, which is equivalent to the integral EZ and the sub- 
sidiary condition (4), defined by 


b 2 
Flu, y, b) = E —-dAK = [ | »( A) + pw? — qyv’ — rv | dx. 


The variation in F resulting from variations in u, y and b is 


>) du (d 6 F ; ; 
5] 2 | y = se) + (p — qy)u du | dx + [ by (4) — qu —dX|dx 
Jo da da J0 da 


du\’ - |* 
+ ly] + (p— qyu —dy| oz # 
dx “ lo 
Inte grating by parts gives 


du (2 iu) du ° [ d ( du) 
ys iJ) a = yu! — du = ly =] dz. 
| “ dx \ dx Y ax lo Jo dx \“ dx 
The integrated part vanishes on account of (6) and é6u = 0 at x = 0. Substituting this 


into the expression for 6F gives 


Pe ee i ( af ae (24 2 
bf —2 | ou! —- y pe — (p — qy)u dx — | by me) = = A | dx 


} 
| 
| a2 / ; 


2 lb 
+ | v(%4) + (p — qy)u’ — wv | éx | . 
: 0 


Of course 6¥ = 0 for any admissible variations du, dy and éx. The two integrands on 
the right give the results already obtained previously. The integrated part vanishes 
at the lower limit because 6x = 0 at x = 0. At the upper limit, 5z ¥ 0. Use of (13) gives 
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y(du/dx)’ + (p — qyw’ — ry = pw =0 at r=b. 
Therefore, 
u=0O at r=b. (14) 
Use of (6), (13) and (14) gives 
y=0 at t=b. (15) 
Because of (2), (1) can be written as 


‘ du dy du (p ™ 
2 = —— = — qy)u. 
"a @eé I ay, 
From (14) and (15), it is seen that 
dy/dx =0 at z= b. (16) 

Thus Eqs. (1), (5), (13), (14) [or (16)] and (15) in conjunction with (4) furnish the 
complete solution for the optimum fin problem. 

CasE (b). Fixed length fin. In practical design, sometimes we wish to restrict the 
length of the fin. If the maximum length permitted is B, the additional restriction is 


b<B. (17) 


The solution from case (a) is valid provided the optimizing value of b satisfies (17). 
Otherwise the solution required is simply for the case of fixed b = B. The procedure 
for affecting a solution is already stated in the first sentence of case (a) except letting 
b = B. The parameter J in u(z, A), y(x, B, A) and H(B, X) is determined from K(B, }). 

3. Heat generation independent of position. If g is assumed to be constant, Eq. (13) 
can readily be solved. For convenience, we define the following notations 


a= ¢@”, 

1+ (1+ r/e)™, 

®(n, —) = Ae *" — de*"/(a’ A), 
®(n, +) = Ae™*” + Ae*"/(a’ A), 
V(n, —) = B(n, —)P(n, +), 

V(n, +) = Are?" + V’e?*"/(a‘ A’). 


ll 


II 


(i) Fixed length fin. From (13) and (5) 
u = 40(z, —) (18) 
use of (18) in (1) with (6) results in 


~B 


y = a [&(2, +)])” | p(x)y¥(x, —) dx. (19) 


z 


The parameter \ is determined from (4) which becomes 


B 
K = }a"(a’ +d)” [ p(x) (x, —) sinh ax dz. 20) 
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Then from (7) 


B 
H = E = halo? +)" / p(av(x, —) de. (21) 
As a > 0, Eqs. (18), (19), (20) and (21) become 

u=1—d'"z, (18a) 

B 
y= [ (\7'”? — x)p(x) dz, (19a) 

B 
K = [ a? - 2)p@r de, (20a) 

aB 
H = | (1 — r'*x)p(x) de, (21a) 


which except for a slight change of notation are the same results obtained from Schmidt's 
criterion [2]. 
(ii) Optimum fin. Equation (18) is still valid. Substituting (14) in (18) gives 


(N/a) optimum = (sinh ab)~’. 


Substituting this relation in Eqs. (18), (19), (20) and (21), and changing B into b give 
the corresponding relations for optimum fins 


u = sinh a(b — x)/sinhab, (22 
b 
y = 3a‘ [cosha(b — x)]~ [ p(x) sinh 2a(b — 2x) dx, (23) 
K = a™’ (cosh ab) [ p(x) sinh a(b — x) sinh ax dz, (24) 
b 
H = (sinh 2ab)™' [ p(x) sinh 2a(b — x) dz. (25) 
It can easily be shown that as a > 0 
u=1—2z/b, (22a) 
y= | (b— 2)p(a) de, (23a) 
K = | (b — x)p(a)x de, (24a) 
H = | (1 — x/b)p(2) de, (25a) 


which are the results from Schmidt's criterion [2]. 

4. Rectangular fin. The variational result just obtained is now applied to specific 
examples of cooling fins. First we consider a fin in the shape of a rectangular plate whose 
sides are of lengths a and b. Let x denote the distance to one of the sides of length a. 
The thickness of the plate is taken to be given by a function y(2). Heat is supplied at 
constant temperature to the side at x = 0. By a suitable choice of units, the temperature 
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at this edge may be taken as unity, and the ambient temperature may be taken as zero. 
Let k be the coefficient of thermal conductivity of the material of the plate, h be the 
combined radiation and convection coefficient of the surface of the plate, and g be 
the constant heat generated per unit volume of the plate per unit temperature change. 
It is then customary to assume that the temperature u(x) satisfies 


d du 
— | ky(x) — | = 2hu — gyu. (26) 
dx | A de v9 
Clearly, Eq. (26) can be written in the form (1) with p = 2h/k = c, where c denotes 
another constant and with g = g/k. Moreover, u satisfies both boundary conditions 


(5) and (6). 

A question in the design of such a fin is the choice of y(x) and b to maximize the 
heat conducted out of the fin base at x = 0 for a given weight fin. On account of (12), 
this heat is 

du 2 


—kay ° =ak | (p — qyudx = akH. 27) 


dx |,-0 Jo 


The volume of the fin is a 
a constant. 
Substituting c for p in Eqs. (23), (24) and (25) gives the solution for the optimum fin 


> y dx = aK. Hence the solution is for the case when p is 


J0 


y = 4ca” tanh’ a(b — 2), (28) 
H = Lea" tanh ab = 3cb ae a K, 29) 
where b is determined from the transcendental equation 
2 


ab — tanhab = 20°K/c. (30) 


If the maximum length permitted is B and if B is smaller than the b determined 
from (30), the solution is then obtained from Eqs. (19), (20) and (21) 


y = hea (W(x, +) — WB, +)]/[&(2, +)), (31) 
H = (1/8)c(a’ + d)7*”[2(2 + A/a’) — W(B, +)], (32) 

where \ is determined from 
(A/d)®(B, +)[e"°* — fala? + dr)’ 8(B, +)] = (4aK/c) — 2B/a. (33) 


5. The circular fin. We now consider a fin in the shape of a circular plate. Let 
t(r) be the plate thickness where r denotes the radial distance. Then the basic differential 
equation is 


d du 
oh Lelie | ee Shes : 
Pe EZ 4 2hru — grtu. (34) 


Heat is supplied to this fin at an inner radius, say r = a. It is then convenient to in- 
troduce the variable x = r — a, which gives the distance to the inner radius. The dif- 
ferential equation becomes 

d 


., du , os _ 
.* | (e + a) | = c(z + alu — g(x + ajlu. (35) 
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The equation is in the form (1) with y = (x + a)t and p = c(x + a). The boundary 
conditions (5) and (6) apply as before if we take the outer radius to be a + b. 
The cooling rate of the fin is given by 2rk f? (p — qy)u dx = 2xkH. The volume 
of the fin is 2x {> y dx = 2K. Again, the solution is given by the relations in Sec. 3. 
Substituting p = c(x + a) in Eqs. (23), (24) and (25) results in the solution for an 
optimum fin 
) = lea? [2(2 + a) tanh’ a(b — x) + a" tanh a(b — x) — (b — 2) sech’ a(b — x)], (36) 
H lea “(1 + 2aa tanh ab — ab sech abe sech ab), (37) 
where b is calculated from 
(b + 2a)(ab — tanh ab) = 4a°K/c. (38) 
If there is a maximum permitted radius, say a + B, then Eqs. (19), (20) and (21) 
can be applied as before. 
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BOOK REVIEWS 


(Continued from p. 238) 


difficulties because standard methods require very small time steps indeed. Alternatively, it is proposed 
to convert the differential equation into a Volterra integral equation and to use numerical quadrature. 

Part IV contains four papers on the numerical solution of partial differential equations. Linear 
parabolic equations are discussed by Herbert B. Keller; the usual finite difference tormulation is given, 
a maximum principle is derived, solv ability established and convergence proved. J. W. Sheldon expounds 
iterative methods for resolving difference schemes arising in the solution of elliptic equations. Accelera- 
tion procedures appropri ite to some of the iterative processes in current use are described and the 
successive overrelaxation method is discussed. Carl N. Klahr next gives a Monte Carlo algorithm for 
elliptic equations first suggested by J. H. Curtiss. There are two papers on hy perbolic equations, one on 
each of the two possible approaches: Mary Lister on that via characteristics, P. Fox on that via finite 
differences. In the first, the mesh is taken in the characteristic, in the second in the original coordinate 
directions—the latter thus requiring careful stability discussions, based on the work of Lax and 
Richtmyer. 

Part V deals with statistical programs. M. A. Elfroysom presents a most useful step-wise multiple 
regression program in which each independent variable is added or deleted at each step, in accordance 
with its contribution to the variance of the dependent variable and with two F values given to the 
program. This procedure will not lead to a unique model, but is insensitive to linear dependence amongst 
the data vectors. Harry H. Harman gives a brief account of factor analysis, with emphasis on a modified 
Jacobi method applied to a reduced correlation matrix. Raymond W. Southworth gives his well-known 
autocorrelation and spectral analysis program (available through SHARE); the autocovariance function 
is computed and raw estimates of the spectral density are found which are then smoothed following the 
ideas of Tukey and Hamming. There is no discussion of aliasing, pre-whitening or trend-removal, nor 
of alternative spectral windows. Lastly, H. O. Hartley presents his method of programming a general 
purpose analysis of variance, which he first described in Biometrika in 1956. 

The final part is entitled ‘‘Miscellaneous Methods” and consists ot the following six papers. Herbert 
S. Wilf treats polynominal equations by a modified Bernoulli method to find a first approximation for 
the root which is then improved by Bairstow’s method. The grave problem engendered by near-multiple 
roots, requiring multiple precision arithmetric, is mentioned and left there. A good discussion of the 
theoretical background and a proof of Sturm’s theorem are included. Anthony Ralston discusses 
Gaussian and Newton-Coates methods for numerical quadrature with a careful consideration of the 
choice of weight-factors; there is no mention of recent work on integration in higher dimensions, but it is 
unfair to criticise individual contributors for omissions since each was clearly not expected to write an 
essay on the current state of his subject but to present a working, and so necessarily limited, computer 
program—and it is this aspect, indeed, which makes the book so exceedingly useful and gives it its 
place as “first of a kind’’. Multiple quadrature by Monte Carlo methods is discussed by Herman Kahn, 
with methods for generating random numbers obeying given probability distributions. In ‘Fourier 
Analysis,’ G. Goertzel evaluates the Fourier coefficients by a simple iterative procedure based on the 
Euler-Maclaurin expansion. This neat method has the advantage that values of only two trigonometric 
functions need be stored or calculated. Dean N. Arden, in a concise but comprehensive and lucid 
exposition of linear programming problems, presents the algorithms of the simplex and of the dual 
simplex methods, with ample theoretical background. The last chapter, by T. R. Bashkow, entitled 
Network Analysis describes the steady-state analysis on a computer of a specific class of ladder networks. 

A book of this sort should be the first of a series: with each article supplemented by a program in a 
recognised algebraic language, such a series would form the backbone of every computer center library. 
The less so can any computer user afford being without the book under review: much programming 
and problem analysis experience has gone into the making of it and this experience is presented here in 
a most palatable and attractive form. 

WALTER FREIBERGER 


(Continued on p. 268) 
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SCATTERING OF A COMPRESSIONAL WAVE BY A PROLATE SPHEROID* 


BY 
NORMAN G. EINSPRUCH anp CARL A. BARLOW, JR. 


Central Research Laboratories, Texas Instruments Inc., Dallas, Texas 


I. Introduction. In order to estimate the size of damaged regions produced in 
single crystal silicon samples by collimated neutron irradiation, Truell [1] utilized a 
multiple scattering of elastic waves calculation made by Waterman [2] to interpret 
ultrasonic attenuation and velocity measurements. The multiple scattering calculation 
assumes a distribution of spherical scatterers and requires detailed knowledge of the 
scattering behavior of a single scatterer embedded in an infinite, homogeneous, isotropic 
medium. The problem of scattering by a single spherical scatterer of a plane compres- 
sional wave which propagates in such a medium has been studied recently by Ying and 
Truell [3]. The analogous problem of scattering of a plane transverse wave has been 
solved by Einspruch, et al. [4]. 

Examination of the single scatterer calculations reveals that the process in which 

portion of the energy in an incident compressional wave becomes transferred into 
energy transported by a scattered shear wave (i.e. mode conversion) is a process of 
higher order; the transfer of energy from an incident compressional wave to a scattered 
compressional wave dominates mode conversion. The equivalent statement for scattering 
of a shear wave holds as well. Since mode conversion can be regarded as a perturbation 
to the scattering process, an approximate description of scattering of a compressional 
wave can be made by replacing the elastic medium in which the scatterer occurs by a 
compressible fluid (i.e. acoustic) medium. 

Since radiation damage in solids [5] frequently occurs as highly localized regions 
which are non-spherical, a prolate spheroid may prove to be a more accurate repre- 
sentation of the shape of a damaged region than a sphere. This is the motivation for 
the solution of the problem presented here. 

In the work of Ying and Truell and of Einspruch, et al. the normalized scattering 
cross section of a spherical obstacle which was a cavity or a mismatched elastic medium 
was found to depend on (Ka)* (K = 2x/\;\ = wavelength of incident wave; a = radius 
of scatterer) in the Rayleigh limit (Ka < 1). In the case of scattering by a rigid fixed 
sphere, the normalized scattering cross section was found to depend on the sum of a 
frequency and size independent term and a term which varies as (Ka)*. Since one would 
expect the scattering cross section to vanish in the limit of large wavelengths, the rigid 
obstacle is probably a poor choice of model for the description of the actual physical 
situation. 

The problem of scattering of a plane compressional wave by a rigid prolate spheroid 
has been solved by Spence and Granger [6]. In the present paper, solutions to the prob- 
lems of scattering by a soft acoustic scatterer and by a small acoustic scatterer of 
arbitrary properties are presented. 

II. Theory. ‘The definition of the prolate spheroidal co-ordinate system is 


x =f cos¢[(# — 1)(1 — 1)]'”, 


*Received December 5, 1960. 
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P(E...) 
(r Op) 











Fic. 1. Prolate spheroidal co-ordinate system. 


y = fsing[(F — IU — 7)]'”, 
2 = jén, 


where f is the semifocal distance, 0 < @ < 27,1 < &< @, —1 < n < 1 (see Fig. 1). 
The mathematical problem is to find solutions of the scalar Helmholtz equation 


V’v+ K*y =0 


within and external to the scatterer which satisfy the boundary conditions at the inter- 
face between the scatterer and the surrounding fluid medium; y is the velocity potential 
function; simple harmonic time dependence is assumed throughout this discussion 

The incident wave is a plane wave of unit amplitude with wave normal in the xz 
plane at an angle © with the positive z axis. The expansion [7] for the velocity potential, 
¥,; , for the incident plane wave is 


oe co . Em : : 
y,=2 p & ys (i)! V.. Sil) Snir (cos O)R'.i(E) cos mg, (1) 
where 
‘) = radial function of the first kind 
St = angular function of the first kind 


Nw = norm of S,,; 
{1, m= 0 
~~ \2,m ¥ 0. 

A thorough discussion of the properties of the prolate spheroidal co-ordinate system 
and the associated wave functions can be found in the text by Morse and Feshbach [7]. 

The radial and angular functions depend on the parameter Kf; wherever it is con- 
venient, this parameter is supressed. The potential function for the wave scattered by 
the obstacle must satisfy the Helmholtz equation, the boundary conditions, and the 
radiation condition 
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pear 
v, ams 1(8, >) r ; 
f(0, ¢) is the scattering amplitude. 
The potential function for the scattered wave is constructed as follows 


He = 2D VO! H™ AmSnu(aRimilE) cos mg, (2) 
wl ere 
A expansion coefficients 
R’ radial functions of the third kind. 


The A,,, are determined by the boundary conditions at the surface of the scatterer. 

III. Scattering by a soft obstacle. Since a soft obstacle will not support a pressure, 
the boundary condition which must hold at the interface between such an obstacle and 
the fluid in which it occurs is 


¥%+t+y¥. = 0, & = 


since the pressure produced by each wave is p 0y/dt, where p is the density of the medium. 
This boundary condition yields 
Rorko 
An: = —Sn: (cos 9) Re mil€o) 


RirilEo). 
The phase shifts are defined as 


(1) 
6.2 = tan” Fe milo) 
—_— < (2) 
R'nilSo) , 
where 
R® = radial functions of the second kind. 
The expansion coefficients are thus given by 


° 4 ~ ibm 
A, = —isin 6,,e '°"'S,. (cos 9). 
Since 
ikKr 
3) /» ° ~\l+m © 
Rig) > — (-)™ Fe 


and at large r 
n — cos 0 
9 d 
erie 2 — ee ee ‘ ae 
{(6, o) —k > yl (—1) yvsin dn; EXP (—75,,2) cos mPS,,; (cos O)S,,; (cos 4). 
‘ ) m=0 + ¥ ml 


The total scattering cross section, o, is defined as 


oa [ [ | {(0, ¢)|? sin 0 dé dé. 
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The integration yields 


8 = = : . 2 -y 7 9 
= ra : pe V sin’ 6,,:[S,: (cos O)}°; 


t=0 m=0 
this expression is exact for all wavelengths. 


In the short wavelength limit 


ae Es l m 
Rni( Kf, &) a Kes” (xe x r) 
and 
es l a 1+ m 
R mi(Af, &) Pp _ Kft cos (if —- 9 = r), 
consequently 


int > — (xe, eS r). 


Kf-2 


The total scattering cross section is thus given by 


Se wes Eo, ee l+m 
oe » ; a [S12 (cos O)]> sin” (x, eimai, 2 


In the long wavelength limit 


m 


Sn. (cos 0) — P" (cos 6) 
Kf—0 


and 


’ 2 (l+~m)! 
N ml = z 
kro 2E+1(1 — m)! 
where the P” (cos 6) are the associated Legendre polynomials. 
The total scattering cross section is thus given by 


+ 1 (J — m)! 


Sn 21 
-_ re p>: p> 2 


IV. Scattering by an obstacle of arbitrary acoustic properties. In the problem 
considered in Sec. III, the common technique of equating the coefficients of the angular 


(1 én SIN’ Sm,[P7 (cos 8)]°. 
+ m)! 


functions which appear in the equation derived from the boundary conditions was used 
to evaluate the expansion coefficients. In the case of the soft scatterer, no wave is excited 
within the scatterer; hence the only waves which propagate are external to the scatterer. 
In the prolate spheroidal coordinate system, the wave number appears as a parameter 
in the angular functions as well as in the radial functions. Consequently, the technique 
utilized previously is not applicable if a wave propagates within the scatterer except 
in the Rayleigh limit, as shall be demonstrated. 

The solution for — > &, is the sum of the incident wave, Eq. (1), and a scattered 


Ss 


wave analogous to Eq. (2). 
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Wi=2>D VGC! 2 Smi(K2, C08 O)Syi(Ko, )Rmi(K2, ) cos mp 


ml 


vy, 2 > pW (i)' 7 A miSmi(Ke, ))Ro(Ke, =) cos me, 


mi 

in these and the following expressions the dependence on the wave number is indicated. 
Subscript 2 indicates a property of the outer medium; subscript 1 indicates a property 
of the scatterer. 

For & < & , the solution is given by 

¥. = 2D DY Gi! Gy BurSmi(Ki, MR ni(Ki, 8) cos mg. 
i=0 m=0 *¥ mi 

In the Rayleigh limit 


vi > 2D LV! y™ PI (cos O)Pi(mRini(K,, 8) cos me (3) 
K,f-0 4 l 


and similarly for y, and y, . 
The boundary conditions which must be satisfied at § = & are continuity of pressure 
and normal component of fluid particle velocity 
po; + pos -_ nel 
= £, e 4 
ay, 4 dy, _— dy. £ = fo 4) 
dé 0& OF ) 


Substitution of Eqs. (3) and Eqs. (4) and performance of the indicated operations yield 





A, = P7 (cos )[pRii(Ko, fo)Rar (Ki, 0) — pind (Ko, &)R(Ki, &)] A"; 
By = P7 (cos O)[pRoi(Ko, fo)Rov (Ke, o) — paRTi(Ke, Ro? (Ko, &)] A”, 
where 
A = pRIMK,, &)RE (Ko, fo) — pa 2(Ka, &o)Ra? (Ky, &0). 
The prime indicates differentiation with respect to &. 
In terms of the expansion coefficients for the scattered wave, the scattering amplitude 
is given by 


{(0,¢) = -F } om >» (—17)” 7 A niSmi(Ke , 9) Cos md. 
mi 


i=0 m=0 4 
As in See. III, the total scattering cross section is obtained by integrating the scattering 
amplitude over a large sphere concentric with the spheroidal scatterer 
re 
m | }2 
o = 72 = | Mee [ - 
K l=0 m + ae 
V. Summary. In this paper the authors have solved the problem of scattering of 
a plane acoustic wave by a soft prolate spheroidal obstacle embedded in a fluid medium. 
An exact solution is given as well as approximate solutions which are valid in the long 
and short wavelength limits. A solution is also presented for the problem of scattering 
of a compressional wave by a small prolate spheroidal obstacle of arbitrary acoustic 
properties. In each case, an expression for the total scattering cross section is derived. 
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ON YIELD CONDITIONS IN GENERALIZED STRESSES* 


BY 
M. SAVE 
Faculté Polytechnique de Mons (Belgium) 


1. Introduction. Generalized stresses and strain rates were introduced in limit 
analysis by Prager [1]. They were defined as the stress-type variables that appear in 
the expression of the power of dissipation, the stress variables being components of the 
stress tensor, stress resultants or dimensionless stress resultants. It was then emphasized 
by Prager [2] that other stress-type variables exist, though not entering in the expression 
of the power of dissipation because they are related to zero strain rates. These stress 
variables he called ‘‘reactions,”’ to distinguish them from those which he called ‘‘gener- 
alized stresses.’’ Numerous applications of these concepts have been made up to now; 
references to most of these can be found in Prager [2] and Hodge [3]. These applications 
essentially concern plate or shell problems, where the generalized stresses are the stress 
resultants and bending or twisting moments. 

The question arises whether it is always possible to express the yield condition in 
terms of stress variables entering in the expression of the power of dissipation, i.e., 
in terms of the so-called generalized stresses only. Physically, the possibility of doing 
this is not obvious because the yield condition essentially involves the stresses, and one 
might think that the non-zero “reactions,’’ because they produce in general non-zero 
stresses, would enter the yield condition. 

We shall show first that it is indeed always possible to express the yield condition 
in terms of the generalized stresses only and, moreover, that the elimination of the 
“reactions’’ can be achieved automatically, by completely ignoring these “reactions.” 

Then, considering the case where one intends to obtain such a yield condition by 
purely statical considerations, we shall prove a theorem that enables one to eliminate 
the reactions from the yield condition in an easy and direct manner. 

Finally, we shall illustrate the theorem by two examples and present some short 
concluding remarks. 

2. Yield conditions in generalized stresses. Let us consider stress variables 
S,, S., +--+ S, used to describe the state of stress of a rigid-perfectly plastic continuum. 


M1, O21 


In general, the yield condition in its normalised form, is 


F(S, , S2, ++: S,) = 1; (1) 
it is represented by a certain hypersurface in the space with the rectangular cartesian 
coordinates S, , S., ++: S,. 

Let us denote 9; , 93, -** g, the corresponding strain rates, and suppose that we have 
qi = 0 ((=kK+1,---n). (2) 

In this ease, S,,, to S, are reactions, while S, to S, are generalized stresses. 

Let the rectangular axes qg; coincide with the axes S; (¢ = 1, --- n). If Drucker’s 


quasi-thermodynamic postulate [4] and its consequences are accepted, conditions (2) 
require that the stress point be on the yield surface (1) at regular points where the 
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normal is perpendicular to all S,., to S, axes or at singular points where this direction 
is on or inside the cone of the outside pointing normals at neighbouring points. 

The points at which conditions (2) are satisfied form a certain hypersurface that, 
projected into the subspace S, , S,, --- S, , gives a yield condition of the form 


GS, , 82, °*°* S&) = 1. (3) 


This projection 1s always possible 

In fact, the only case of impossibility would occur if at least one point determined 
by (2) on the yield surface (1) would be on one of the S,,, to S, axes, say S,; . The pro- 
jection which, strictly speaking would still be possible, would give, instead of a certain 
set of S, , S.,--- S, satisfying (3), the value S, = S, = --- = S, = 0. 

But this case is impossible because, due to the convexity of the yield surface, this 
would put the origin on or outside the yield surface, which is not admissible. 
the stress variables entering 


This being shown, we can call ‘‘generalized stresses’ 
in the expression of the power of dissipation or in the yield condition, denote them 
by Q, ,Q., --- Q, and re-write (3) as 


(Q, ,Q., °° Q,) = 1. (3’) 


The yield condition (3) can be obtained directly, ignoring completely the “‘reactions,”’ 


as it was done in [5] and elsewhere. 
To show this in general terms, let us return to (1) and consider the power of dissipation 


D= Si: Gi “ So* Qs + eee + Sid: F (4) 
The perfectly plastic continuum being inviscid, D is an homogeneous function of 
order one of q; , -:- g, and we can write 
aD , oD . 
D — = 7h ~_ “ee — = A) a i ep? 
07”: O”d, 
Comparing (4) to (5) we obtain 
aD 
; == (2 = 1,2, ---n). (6) 
0”; 


Equations (6) are a parametric definition of (1), the parameters being the g; (or 
their ratios). 
Now if we have conditions (2), (4) reduces to 
D*¥=Q:-qGi + QeGte Hau, 
the stress variables S, to S, being generalized stresses which we denote consequently 
Q, to Q, . The yield condition (3’) is given by the first k equations (6) where conditions (2) 
are entered. 
But, as only partial differentiation is used, we obtain the same relation in the q; 
(i = 1, --- k) if we form dD/dq,; (¢ = 1, --- k) and use g O(@ =k+1,--+n) or 
if we form dD*/dq, 
Thus, the first k equations (6) are identical to 
aD* 
09; 


(ta 1, --- k). (7) 
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Equations (7) give directly (3’) under parametric form, ignoring completely all 
reactions S, 
We shall now discuss the following two points: 
|. Does the normality law apply to the yield condition (3’) and, more generally, 


to S, 


l 


in arbitrary linear subspaces? 

2. Is it still correct to ignore ‘‘reactions’’ when the yield condition is established 
by purely statical considerations or, in other words, how can reactions be elimi- 
nated from (1)? 

3. Normality law in subspaces. The validity of the normality law in general 
subspaces is by no means obvious and should be investigated first. 

We denote by S, , --- S, all non-vanishing stress variables. The yield condition 
in its normalised form is (1) where all non-vanishing stress variables enter, in general. 

In the most general way, one obtains a certain simplified yield condition in a certain 
subspace when there exist a certain number of relations between the stresses, say (n — k) 


relations, 


f,(S,,---S,) =0 

fo GS, ’ S, on 0 

by : (8) 
2 

lf, AS, silt S,) = (0). 


If we suppose that we are able to express, by means of (8), (n — k) stresses as func- 


tions of the others, Eqs. (8) can be re-written 


|S, 1 = SpiilS, , &, S,) 
’ ’ ’ ’ 
1Ses2 = Ses2(Si , Se, S,) (9) 
S, Be ie 5 Bes S,) 
Introduction of (9) in (1) furnishes a simplified yield condition 
WS, ates S,) = | (10) 


in the subspace (S, , --- S,). 
We consider first regular yield surfaces. 
The normality law, applied to the initial yield condition (1) gives the strain rates 


AF . 
qi =< (i = 1,2, ---n), (11) 


\ being a positive constant. 
The normality law remains valid for the yield condition (10) in the subspace 


(S, , --- S,) if we can write* 
ov 

ad (= 1,-+-®. 12 

q Os, ~~ 


*q: fori = | 1, --* n cannot be obtained of course. 
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We therefore have 








oF ov : 
as, _ as, (a = a; pli k). (13) 
On the other hand, we know that 
W(S, , S&S) = F[S, pes S;, ’ SiiGS, anita Sx); sc SCS, aes S,)] (14) 
Applying chain differentiation, we find 
ov OF OF OS,.; OF OS,, ° pers 
88, 9S, dSia OS, * °° +98. as, “=i (19) 


In general we cannot expect relations (8), or their explicit form (9), to be such that 
(15) reduces to (13) and so ‘‘in general, the normality law is not applicable to simplified 
yield conditions.” 

This law will however continue to be valid in two important cases: 

a) If relations (9) have the special form 


S; = K; (2 = k + ie a n), (16) 
where the K; are constants (such, of course, that the yield condition is not 


violated). 
In this case we have 


ee OS, +2 OS, . yj 
SS es ee ee ©, (17) 
ON; ON; OS; 


and (15) reduces to (13). 
b) If S,41 to S, are “reactions.” 
In this case we have 





q; = 0 (@=k+1,---n) (18) 
and, applying the normality law (11), relations (8) become 
oF ; 
—— = = » “ee 
as. = 0 @=k+1, n) (19) 
or, explicitly, 
a a (20) 


a. teas oe 


and (15) reduces once again to (13). 

We may sum up our findings as follows. 
“The normality law remains applicable to simplified yield conditions obtained by giving 
constant values to a certain number of stress variables or by eliminating the reactions.” 

Let us now turn to yield surfaces with singular points. The normality law is gener- 
alized by the condition that the strain rate vector (q;) be in the cone of the outside 
pointing normals at the neighbouring points of the pointed vertex. 

Equation (10) represents the cylinder with generating lines parallel to the axes of 
Sis: to S, , projecting the intersection of (1) with the hypersurfaces (8). In the space 
S,, °*: S,, it is the equation of the cross section of this cylinder. 
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The normality law will remain valid for the hypersurface (10) if, at a pointed vertex 
P of (1) that is also on (8), the projection of the cone of normals at (1) in P is the cone 
of normals at (10) at the point P’ projection of P. 

We consider now a pointed vertex as the limit of a very small region where the 
direction of the normal varies very quickly but continuously. 

Suppose that this region is bounded on the yield surface by a certain curve C, whose 
limit is the point P at the vertex. When C tends to P, we can at each stage express 
Eq. (13) in the following terms: the projection of the cone of normals along C are the 
normals to the projection of C. 

Considering the conditions of reduction of (15) to (13) and passing to the limit, we 
readily extend our preceding conclusions to singular yield surfaces. 

4. Theorem on the adaptation of the reactions. Let us now suppose that, among 
the totality of nm non-vanishing stress variables, the k first ones are generalized stresses 





Q, ,Q., +: Q, and the others are reactions S,4; , Sis2, -** Sn- 
We give fixed values to all generalized stresses but one, say Q, 
Q; = K, 
ata (21) 
Q.-1 aa Ky-1. 
We obtain a simplified yield condition 
9(Q, , Sear, °** S,) = F(K, , Ko, ++ Ku-s » Qe , Sear, -°* S,) = 1 (22) 


when we use (21) in (1). 
Suppose that the yield surface is regular. 
The normality law being valid for the yield condition (22), we have 


fefte) 





oe. _ scoiall 9 
Gi=A5g = RAI, +n) (23) 
(and also g; = 0¢/0Q, which is of no interest for our purpose). 
The stress variables S,,, to S, being reactions, 
qi; = 0 (@=k+1,--+n). (24) 
Comparing (23) and (24) we obtain, 
00 ° . ~ 
— = 0 (2 =k +. 1, coe n), (25) 
OS; 


because A + 0. 

Equations (25) are extremum conditions for Q, , considered as a function of 
s+. , *** S, implicitly defined by (22). 

Let us note, that, as the yield surface is convex at all points, relations (25) are con- 
ditions of absolute maximum or minimum. 

We can then state the adaptation theorem: 

“Tf we fix all generalized stresses but one, the reactions adapt themselves in such a way 
as to give to the non-fixed generalized stress a maximum positive value or a minimum 
negative value.” 


S 
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The theorem remains valid for a singular yield surface as can be readily shown: 
the generalized normality law applies to the yield surface (22), and conditions (24) 
determine those points on (22) where the strain rate vector is parallel to the Q, axis. 
If this occurs at a pointed vertex, Q, assumes an extreme value at this point because 
of the convexity of the yield surface, which contains the origin. 

5. Examples of application. We consider a rigid-perfectly plastic cylindrical shell 
without axial load, under axially symmetrical loading. As in most plate and shell prob- 
lems, we suppose that straight lines normal to the undeformed middle surface remain 
normal to the deformed surface during yielding. The axis 0Z directed along the normal 
at a typical point of the middle surface of the shell (Fig. 1) is then a principal axis of 





the strain rate tensor. Consequently, as Ziegler [6] has shown, in the case of the Tresca- 
Guest yield condition as well as in the case of the von Mises yield condition, it is also 
a principal axis of the stress tensor and 7,, = 7., = 0 at all levels in the thickness of 
the shell. Consequently, to be consistent with the previous hypothesis, shearing forces 
normal to the middle surface in a completely plastic element of the shell (or of the plate) 
must vanish exactly. Moreover, we suppose ¢, = 0. We are then exactly in plane stress 
conditions. 

Due to the symmetry of our problem, the only non-vanishing stress variables are 
M, , M, and N, . (Fig. 2) 








Fic. 2 


Symmetry also imposes the condition that the circumferential rate of curvature 
Xe be zero. 

M, being then a reaction, we have to determine the yield condition in M, and N, . 

1. Let us first consider a sandwich shell, obeying von Mises’ yield condition. This 
shell is composed of two thin sheets, each of the thickness t/2, which are separated by 
a core of the thickness h. (Fig. 3). The thin sheets support only normal stresses, uni- 


formly distributed on the thickness. 
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If we eall 
o,, and o,, the stresses in the top sheet 
o», and a,, the stresses in the bottom sheet 





we have 
( t 
= 5 (on + om), 
lat, = © (on — om) (26) 
| 
M, = f(a, — G2). 


N, _M. M, 


No = —; Ms , 
M, 


where N o,-tand M, = o:th/2, 


o, being the yield stress of the sheets in pure tension, 





we get 
Oar + Onan 
te = ere 
<0 
o — Cap, Qo” 
yo (27) 
205 
G24 —~ Gy 
mm, =~ ao ° 
“04 
von Mises’ yield condition 
ao, + 0 — 0,09 = Go (28) 


is shown in Fig. 4 in the (Oo, , 0o,) plane. In this plane, the state of stress in the shell 
is represented by a straight segment ¢b, coordinates of ¢t and b being respectively (¢, , z+) 
and (a9, , ¢ In view of (27), coordinates of the center c of tb give ny (and n,) and 
projection of the segment tb on the axes give m, and m, (the positive factors 1/0) and 


Jf 


1 /2c, do not create any trouble). 

As n, 0, the point ¢ must stay on the Oo» axis. For a given position of this point 
(between A and B), corresponding to a given value of ng , plasticity of the shell element 
requires one point ¢ or b at least to be on the yield locus. 

The adaptation theorem tells us now that the inclination of t must be taken so 
that its projection on Oc, be a maximum. This readily gives the condition 
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m 
Mm = > (29) 


~ 


which moreover puts both points, ¢ and b, on the yield locus. Displacing the point c¢ 
from A to B and respecting condition (29), we obtain without difficulty and directly, 
the elliptic relation 


3m: +n, =1 (30) 


between m, and n, . 
2. Let us consider the same shell, but obeying now Tresca-Guest’s yield condition. 
A very similar analysis, with the use of the adaptation theorem gives the following 
results (Fig. 5) 


a. when 0 < n < 3, 





m 
Me <—<1—m and m, =1 (31) 
mM, 
b. when 3 < my < 1, ms = m,/2 and 
Io 
oA b 
t 
ed Oz 
Vo weir ace 
15 


Fic. 5 
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= +n =1 (32) 
2 

The yield condition is given by (31) and (32). 

6. Concluding remarks. As has been shown in Sec. 3, the adaptation theorem is 
by no means the only way to obtain the correct yield condition in generalized stresses. 
Even in the case of a purely statical procedure, reactions can also be eliminated by 
other considerations. 

It should be pointed out, however, that the maximizing procedure applied in Sec. 5 
has already been used (see [7] as an example). 

The method was open to criticism at that time because it was not certain that re- 
actions, which were ignored, would in fact take values which permitted this maximizing. 

Our theorem gives, a posteriori, a rigorous justification of this method. 

Acknowledgment. The author is indebted to professor W. Prager for valuable 
suggestions. 
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Ludwig Prandtl. Gesammelte Abhandlungen zur angewandten Mechanik, hydro- und 
Aerodynamik. Edited by W. Tollmien, H. Schlichting, H. Gértler, and F. W. Riegels. 
Vol. I, xix + 574 pp., Vol. I, xiii + 495 pp., Vol. III, xiii + 550 pp. Springer- 
Verlag, Berlin, 1961. $68.96 (paperbound), $72.75 (clothbound). 


This publication of the collected works of Ludwig Prandtl has been sponsored by the Max-Planck- 
Gesellschaft zur Férderung der Wissenschaften, the Wissenchaftliche Gesellschaft fiir Luftfahrt, and 
the Gesellschaft fiir Angewandte Mathematik und Mechanik. In addition to an excellent portrait 
of Prandtl, a Preface by A. Betz and Th. v. Kaérm4n, an Introduction by W. Tollmien, and a chronological 
list of publications, the first volume contains contributions to elasticity, plasticity, and rheology, as well 
as wing and propeller theory. The second volume is devoted to papers on boundary layers and drag, 
turbulence, creation of vorticity, and gasdynamics. The third volume contains papers on meteorological 
applications, model experiments, and miscellaneous subjects, as well as lists of publications that have not 
been reproduced in the present work, of dissertations written with Prandtl’s guidance, and of bio- 


graphical data and awards. 


Computing methods and the phase problem in X-ray crystal analysis. Edited by Ray 
Pepinsky, J. M. Robertson and J. C. Speakman. Pergamon Press, New York, Oxford, 
London, Paris, 1961. viii + 326 pp. $9.00. 


This is a collection of 28 papers and ensuing discussion presented at a conference held at Glasgow 
in August, 1960. A similar volume under the same title was issued in 1952 by the first-named editor; 
the present work is, in a sense, volume IJ. The international group of authors of the papers include 
a large fraction of those active in more extensive crystallographic computations and in general attacks 
on the phase problem. 

Somewhat more than half of the book is devoted to the strategy of crystallographic calculations 
on a variety of computers. The calculations include evaluation of three-dimensional Fourier summations, 
calculation of structure factors, and refinement of trial structures by least-squares and other methods. 
Computer programs are described in general terms so that the methods may be widely applied. This, 
plus description of much practical experience, should make this book extremely valuable to someone 
who wishes to program similar calculations on a new computer. 

What is called the ‘‘phase problem”’ in crystal analysis arises from the fact that the experimental 
method gives the magnitudes, but not the phases, of Fourier coefficients needed to work out the structure 
of a crystal. Methods of solution range from trial-and-error to highly sophisticated statistical approaches. 
This volume constitutes a progress report on new developments in this field, which seem to be less 
numerous and less optimistic than those of five years ago. 

This up-to-date and authoritative book is recommended to crystallographers and to computer 


experts who deal with them. 
G. B. CARPENTER 


Elements of mathematical statistics. By Howard W. Alexander. John Wiley & Sons, 
Inc., New York, London, 1961. xi + 367 pp. $7.95. 


This book is intended as an introduction to probability and mathematical statistics for students 
who have completed a year of calculus. A feature of the presentation is that even at this elementary 
level probability theory is erected on a foundation of set theory. This gives the student the proper 
orientation for more advanced work. 

WALTER FREIBERGER 
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—NOTES— 


A NOTE ON AN INVERSE PROBLEM IN MATHEMATICAL PHYSICS* 


BY 
RICHARD BELLMAN (The RAND Corporation) AND 
JOHN M. RICHARDSON (Hughes Research Laboratories) 


1. Introduction. The problem of solving the linear partial differential equation of 
diffusion theory and heat conduction, 


7 2 
du _ ou (1) 


subject to appropriate initial and boundary conditions has been extensively investigated 
using a variety of analytic techniques. We wish to consider a type of problem associated 
with this equation which has not been investigated to any extent, the inverse problem. 

In place of assuming that the coefficient function g(x) is known and trying to find 
the solution u(x, t), we shall suppose that the solution is known at certain points of the 
x, t plane, and try to determine the function ¢(x) from this information. This type of 
question appears first to have been asked by the mathematical physicist Ambarzumian, 
and it was also posed by Langer, who gave an analytic procedure for determining the 
coefficient function. The first detailed analysis of this problem was given by Borg [1], 
where references to earlier work by Langer and others may be found. Similar questions 
arise in the quantum theory of scattering, and were, in principle, solved by Gelfand 
and Levitan [2]. The most recent paper is that of Kay [3], where references to work of 
Levinson and others may be found. 

The aim of this paper is to show how problems of this nature are closely related to 
some classical synthesis problems for electrical networks. Our new approach furnishes 
a way of obtaining approximate solutions which appears to be more useful from the 
computational, engineering and physical point of view than the exact solutions furnished 
implicitly by the methods of Borg and Gelfand and Levitan. 

2. The analytic problem. Consider the equation of (1.1) taken to hold for0 <a < L, 
t > 0, subject to the initial and boundary conditions 


(a) u(L, t) = 0, 
(b) u(x, 0) = 0, (1) 
Ou 





az leno > | 2(4 uc |. ; 


Here Z(s) is a meromorphic function of s which we shall specify in further detail below. 
Taking Laplace transforms with respect to t, we obtain the Sturm-Liouville problem 


(a) initia: se ee 

\a oY ye _ dx?’ 

(b) o(L, s) = 0, (2) 
dv ‘ 

ad dx lr=0 oo Z(s)u r=0 : 
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where 


--) 


v(x, 8) = | u(x, te ** dt. (3) 


J0 


Given the function Z(s), we wish to determine the density function. For an ap- 
propriate choice of Z(s), this problem is equivalent to the problem considered by Borg. 
To see this, suppose that Z(s) has the form 


Zs) = [Ja +s8/,,)/J] 0 + s/n), (4) 
n=1 n=1 
where X,, , #, > 0. If s = —X, , we have the Sturm-Liouville problem 
v’ + rv = 0, 
v’(0) = o(L) = 0. 
If s = —y, , we have the problem 


v’ + uw = 0, 
v(0) = o(L) = 0. 


(6) 


We thus have the two distinct problems required for the solution of our problem 
following Borg’s methods. Unfortunately, Borg’s results do not furnish any useful 
technique for actually obtaining the unknown function, g(x). We shall present another 
approach which thus furnishes an alternative attack upon Borg’s problem. 

3. Discrete approximation. Let us, following a standard route, replace the differential 
equation of (2.2a) by the difference equation 


Vrs — 2, + Yy:-1 = A’sy,r, , r= Z.- es N= §, (1) 
where 
v, = v(k A), o& = ofk A), BN A e-L. (2) 
The boundary conditions of (2.2b) and (2.2c) are replaced by 
a vy = O, e 
(a) N (3) 
(b) Vv, — Uo = AZx(s8)v . 
We replace the meromorphic function Z(s) by its approximation 
II a +s) 
Zx(s) = Treen . 
IT (1 + 8/n,) 
n=1 
4. Network theory. Once the problem has been transformed into that of deter- 
mining the values of g, , an approximation to the determination of g(x), we see that it 
is equivalent to that of determining the elements of a two-port network, a network of 
transmission-line type, given the response function. 


Probably the most efficient and convenient way of deriving the values of the ¢, 
is based upon the use of continued fraction expansions. The details of this procedure, 


(4) 
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together with various necessary and sufficient conditions upon Z(s) and the zeros and 
poles, —\,, and —u, , may be found in Guillemin [4], Weinberg [5], Weinberg and Slepian 
[6], where many further references may be found. 

Consider the system of equations 


v 1 — (2 — A’sy,)y, + %-1 = 0, e@ 1,2, +++, 87 — 1, () 
Writing 
v;- 42 / / 
“th, = (2 ~- & Spr) el 1, (v/Ve+1), (2) 
k 


we see that we obtain the finite continued fraction 





. 1 
= (2 — A’sy,) — 2 
V; , s¢1) (2 — A’sp2) — ; 


(3) 


Hence, any method which permits us to write vo/v, , as given by (3.3b), as a continued 
fraction of this form, yields the values ¢, , g2 , «++ , ¢v , and thus gives us a computational 
solution of the discrete inverse problem. 

Questions of degree of approximation will be discussed subsequently. 

That there is a strong connection between the mathematical problems of network 
theory and those of quantum mechanics is no longer surprising in view of the work of 
Wigner and others on positive real functions. See Wigner [7], Lane and Thomas [8]. 
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